package functions;

import svc.BioRecordSvc;
import svc.ResultSvc;
import utils.Utility;
import constants.Constants;
import domain.BioRecord;
import domain.DelRecord;
import domain.EstRecord;
import domain.FutRecord;
import domain.HRecord;
import domain.InsRecord;
import domain.MkvRecord;
import domain.PolRecord;
import domain.RpoRecord;
import domain.TestRecord;

public class Functions {
	BioRecordSvc bioSvc = new BioRecordSvc();
	ResultSvc resSvc = new ResultSvc();
	
	public void startSimulation(int fileNo){
		BioRecord bio = new BioRecord();
		FutRecord fut = new FutRecord();
		PolRecord pol = new PolRecord();
		InsRecord ins = new InsRecord();
		RpoRecord rpo = new RpoRecord();
		DelRecord del = new DelRecord();
		MkvRecord mkv = new MkvRecord();
		
		bio = bioSvc.getBioRecordDetails(fileNo);
		
		bio = bio_probability(bio);
		int x1 = dp_policy(bio, fut, pol, ins, rpo, del, mkv);
		System.out.println("x: "+x1);
		
		int n=0;
		double conv = 0.0, conv1 = 0.0;
		
		mkv.setStart_day(1);
		mkv.setStart_sea(bio.getStart_sea());
		mkv.setStart_cows(0);
		
		
		if(bio.getAnalysis() == 1)
		{
			do													//calculate steady-state herd
			{
				if(n > 1) conv1 = conv;
				for(int i = 0; i <= Constants.CURR; i++)						//curr
					for(int j = 0; j <= bio.getMaxlact(); j++)			//lact
						for(int k = 0; k <= Constants.MIM; k++)				//mim
							for(int m = 0; m <= Constants.MPREG; m++)			//mpreg
								mkv.getHlp()[i][j][k][m] = mkv.getCow()[i][j][k][m][2];
				//for(i=1;i<=5;i++)    //this would have error calcutaled every 5x12=60 months, but appears to be hardly faster
				markov_chain(1, Constants.SEA, 0, bio, mkv, pol, rpo, ins, del);							//1=DP, SEA=12 months, 0=do not calculate statistics
				conv = 0.0;										//take sample of new markov chain
				for(int i = 0; i <= Constants.CURR; i++)						//curr
					for(int j = 0; j <= bio.getMaxlact(); j++)			//lact
						for(int k = 0; k <= Constants.MIM; k++)				//mim
							for(int m = 0; m <= Constants.MPREG; m++)			//mpreg
								conv += Math.abs(mkv.getHlp()[i][j][k][m] - mkv.getCow()[i][j][k][m][2]);
				System.out.println("   steady-state conv:"+conv+" ("+(++n)+")");
				
				if(conv > conv1 && n > 100) 
				{ 
					System.out.println("convergence decreases, continue with current solution"); 
					conv = 0.0;
				}
				if(n == 500) 
				{ 
					System.out.println("steady state not found in 500 iterations, continue with current solution"); 
					conv = 0.0;
				}
			}
			while (conv > 0.0000001);				//convergence: herd in steady state

			markov_chain(1, 36, 1, bio, mkv, pol, rpo, ins, del);					//calculate statistics for 4 years
			markov_chain_stats(36, bio, mkv);					//calculate markov chain statistics
		}
		
		resSvc.saveAllOutputs(fileNo, bio, mkv, pol, rpo, ins, del);
		//WriteFunction.writeFunction(36, bio, mkv, rpo, pol);
		
	}
	
	@SuppressWarnings("unused")
	public BioRecord bio_probability(BioRecord bio)
	{
		int sea = 0, curr = 0, next = 0, lact, mim = 0, ok = 0, run = 0, lastrun = 0, i = 0;
		double  r21		= 0,	//0.303,	//correlation coefficient with previous lactation weight (0.55^2)
				rc21	= 0,	//.550,	//regression coefficent next on curr (1st lactation)
				sd		= 0.0,		//coefficient of variation of prediction
				sumprob = 0.0,
				beststep= 0.0,
				diff	= 0.0,
				maxdiff	= 0.0,
				step	= 0.0,
				goal	= 0.0,
				m305	= 0.0,
				xlm = 0.0, xum = 0.0, zxlm = 0.0, zxum = 0.0, pxlm = 0.0, pxum = 0.0;
		double b21[] = new double[Constants.MIM+1];
		double px[] = new double[Constants.CURR+1];
		double funk[] = new double[Constants.MIM+1];
		double yum[] = {0, 74,78,82,86,90,94,98,102,106,110,114,118,122,126,130}; //upper limits of levels
				//Funk adjustment factors for previous days open on current 305-d milk yield. See Funk[1200] and
				//   valueofpregnancy.xls!funk
	
		System.out.println("\n BIO_PROBABILITY started ...");
		if(bio.getPa_cv() == 0)				//special case: no transition to other classes
		{
			bio.getPnext()[0][0][0][0][8] = 1.0f;
			for(lact = 0; lact <= Constants.LACT; lact++)
				for(mim = 0; mim <= Constants.MIM; mim++)
					for(sea = 0; sea <= Constants.SEA; sea++)		
						for(curr = 0; curr <= Constants.CURR; curr++)
						{
							bio.getPnext()[curr][lact][mim][sea][curr] = 1.0f;
							bio.getLo()[curr][lact][mim][sea] = curr;
							bio.getHi()[curr][lact][mim][sea] = curr;
						}
			return bio;
		}
	
	
		//if b21 == 1 when mim >= 1 and b21 = 0.55 when mim == 0 then we are back to the
		//old model where cows only move from lactations.
		//bio.getpnext[0][0][0][next]        : probability heifer is in a 'next' level in lactation 1
		//bio.getpnext[curr][lact][mim][next]: probability cow is in 'next' level in next month
		//b21[] correlations taken from Houben[816]2986.
		//
		//August 8, 2005:  MILK TRANSITION PROBABILIES and EFFECT OF PREVIOUS CALVING INTERVALS
		//Houben[816]2986 uses both current and previous lactation. 
		//It is not clear how Houben uses the method of Van Arendonk[128]118 with the correlations in Houben[816]2986
		//He calculates this in DPAMELK.PAS, but this is too complex to understand.
		//Groenendaal[]2150 uses 0.55 and 0.50 but they have only the current lactation. 
		//    Maybe they predict milk production in the future instead of backwards like in DP? 
		//Van Arendonk[128]119 uses 0.55 and 0.50 but he has both current and previous lactations.
		//Van Arendonk[129]16 uses 0.60 when he uses only the current lactation.
		//I use only one lactation, but still want to use milk level transition probabilities between months.
		//It is not clear how to do that from Houben's paper.
		//Let's assume that Constants.CURR level is cumulative level of milk production from mim=1 to mim = now.  Then when mim progresses,
		//   we are getting more sure about the current level of milk production.  So Constants.CURR is the predicted level of 305 d milk?
		//-.get I am going to use 0.60 according to Van Arendonk[129]16 
		//Van Arendonk[674] + Dekkers[490]239 include effects of previous calving interval on current milk production.
		//   I would like to add this: can be done when cow moves to next lactation.
		//Effect of number of days open in current lactation is included through effect of pregnancy on milk yield (
		//I would also like to add the effect of previous lenght of lactations
		//-.get shift distribution of transition to next lactation, use factors in Funk[1200] 
	
		for(mim = 0; mim <= 10; mim++) funk[mim] = 1.0;
		funk[11] = 0.980392157;	//2 months open -.get calving when mim == 11, Constants.MPREG == 9
		funk[12] = 1.000000000; //3 monhts open -.get calving when mim == 12, Constants.MPREG == 9
		funk[13] = 1.010101010;	//Funk[1200]
		funk[14] = 1.020408163;
		funk[15] = 1.030927835;
		funk[16] = 1.036269430;
		funk[17] = 1.041666667;
		funk[18] = 1.047120419;
		funk[19] = 1.052631579;
		funk[20] = 1.052631579;
		funk[21] = 1.052631579;
		funk[22] = 1.052631579;
		funk[23] = 1.052631579;
		funk[24] = 1.052631579;
	
	
		b21[0] = bio.getMilkrepeat();
		b21[1] = 0.956;	//correlation between total milk yield in month 1 and total milk yield in month 2
		b21[2] = 0.979;	//correlation between total milk yield in month 2 and total milk yield in month 3
		b21[3] = 0.988; //month 3 and month 4
		b21[4] = 0.991; //month 4 and month 5
		b21[5] = 0.993; 
		b21[6] = 0.994;
		b21[7] = 0.996;
		b21[8] = 0.997;
		b21[9] = 0.997;
		for(mim = 10; mim <= Constants.MIM; mim++) b21[mim] = 1.0;	//cow remains in same curr level
	
		if(bio.getFunk() == 1) bio = predict_305m(bio);				//predict 305-d milk production
		
		for(lact = 0; lact <= Constants.LACT; lact++)
		{
			for(mim = 0; mim <= Constants.MIM; mim++)
			{
				for(sea = 0; sea <= Constants.SEA; sea++)				//sea == 0: transition within lactation
				{											//sea >= 1: transition between lactations
	
					if(lact == 0)							//heifer classes, no previous information
					{
						rc21 = 0.0;
						r21 = 0.0;
					}
					else if(sea >= 1)						//cow classes: transition to next lactation
					{
						rc21 = b21[0];						//regression coefficient = repeatability
						r21 = b21[0] * b21[0];				//correlation coefficient = repeatability^2
					}
					else	//sea == 0						//cow classes: transition within lactation
					{
						rc21 = b21[mim];					//regression coefficient = repeatability
						r21 = b21[mim] * b21[mim];			//correlation coefficient = repeatability^2
					}										//b21[0] is transition to next lactation
					sd = bio.getPa_cv() * Math.sqrt(1 - r21);
					for(curr = 1; curr <= Constants.CURR; curr++)		//calculate P(next level | current level)
					{
						ok = 0;
						run = 0;
						step = 0.0;
						lastrun = 0;
						while(ok == 0)
						{
							m305 = 0.0;
							for(next = 1; next <= Constants.CURR; next++)
							{
								if((rc21 < 0.9999 || lact == 0)) 
							//	&& bio.getSmy()[curr] > 0)	//if no bio.getSmy() constraint then weird results for small sd
								{		//problems would occur when b21[mim] ~= 1.0
									if(next == 1) xlm = -9999; else xlm = xum;	//z-score
									xum = (yum[next] - rc21 * (bio.getSmy()[curr]-100) - 100)/sd;
	
	
									if(next == Constants.CURR) xum = 9999;
									zxlm = 1/Math.sqrt(2*Constants.PI)*Math.exp(-0.5*Math.pow(xlm,2));			//not needed: P(z-score)
									zxum = 1/Math.sqrt(2*Constants.PI)*Math.exp(-0.5*Math.pow(xum,2));			//not needed: P(z-score)
									pxlm = Utility.normal(xlm+step);							//cumulative P(z-score)
									pxum = Utility.normal(xum+step);
									bio.getPnext()[curr][lact][mim][sea][next] = pxum - pxlm;	//P(level 'next' | level 'curr')
	
	
								}
								else	//b21[mim] ~= 1.0
								{
									if(next == curr)
										bio.getPnext()[curr][lact][mim][sea][next] = 1.0;//P(level 'next' | level 'curr')
									else
										bio.getPnext()[curr][lact][mim][sea][next] = 0.0;
								}
	
	
								if(lact >= 1 && lact <= 3 && mim >= 11 && mim <= Constants.MIM && sea >= 1)	//transition to next lactation
								{
									m305 += bio.getPnext()[curr][lact][mim][sea][next] * bio.getM305()[next][lact+1][sea];
							
									if(lastrun == 1 && lact == 3)
										for(i = 4; i <= Constants.LACT; i++)
											bio.getPnext()[curr][i][mim][sea][next] = bio.getPnext()[curr][i][mim][sea][next];
								}
							}//next
	
							if(bio.getFunk() == 0) ok = 1;
							else if(lastrun == 1) ok = 1;
							else if(lact >= 1 && lact < 3 && mim >= 11 && mim <= 24 && sea >= 1)
							{
								if(run == 0) 
								{
									goal = m305 * funk[mim];
									maxdiff = 100000;
									beststep = 0.0;
									step = -2.0;
									px[0] = m305;
									for(next = 1; next <= 15; next++) px[next] = bio.getPnext()[curr][lact][mim][sea][next];
								}
								else
								{
									diff = Math.abs(m305 - goal);
									if(diff < maxdiff) beststep = step;
									if(diff < maxdiff) step = step + 0.05;
									if(diff >= maxdiff || run == 80)	//no improvent -.get run once more with beststep;
									{
										step = beststep;
										lastrun = 1;
									}
									if(diff < maxdiff) maxdiff = diff;
								}

							}
							else ok = 1;
	
					
							run++;
						};//while
	
	
						sumprob = 0.0;
						for(next = 1; next <= Constants.CURR; next++)
						{
							if(bio.getPnext()[curr][lact][mim][sea][next] < Constants.PRUNE) 
								sumprob += bio.getPnext()[curr][lact][mim][sea][next];	
							else
							{
								bio.getLo()[curr][lact][mim][sea] = next;
								break;
							}
						}
						for(next = Constants.CURR; next >= 1; next--)
						{
							if(bio.getPnext()[curr][lact][mim][sea][next] < Constants.PRUNE) 
								sumprob += bio.getPnext()[curr][lact][mim][sea][next];	
							else
							{
								bio.getHi()[curr][lact][mim][sea] = next;
								break;
							}
						}
		
	
						for(next = 1; next <= Constants.CURR; next++)
						{
							if(next < bio.getLo()[curr][lact][mim][sea]) 
								bio.getPnext()[curr][lact][mim][sea][next] = 0.0;
							else if(next <=  bio.getHi()[curr][lact][mim][sea])
								bio.getPnext()[curr][lact][mim][sea][next] = bio.getPnext()[curr][lact][mim][sea][next] / (1 - sumprob);
							else
								bio.getPnext()[curr][lact][mim][sea][next] = 0.0;
					
	
							//bio.getpnext[0][0][0][0][next] are the heifer transition probabilities but curr == 0 would not be 
							//filled after moving the curr classes from 0-14 to 1-15.  
							if(lact == 0 && mim == 0 && curr == 8)
								bio.getPnext()[0][0][0][0][next] = bio.getPnext()[8][0][0][0][next];
						}//next
	
						sumprob = 0.0;							//check: do probabilities sum up to 1
						for(next = 1; next <= Constants.CURR; next++) 
							sumprob += bio.getPnext()[curr][lact][mim][sea][next];
						if(sumprob > 1.000000000001 || sumprob < 0.999999999999) 
						{ 
							System.out.println("\n BIO_PROBABILITY sumprob:"+sumprob+" != 1.000"); 
							Utility.stop(1);
						}
	

					}//curr
				}//sea
				if((mim == 8 || mim == 16) && lact <= 3 && lact >= 1)
				{
					System.out.println((double)(lact-1)*6+(double)mim/3);
					Utility.write_msg_file(22,(double)(lact-1)*6+(double)mim/3);
				}
			}//mim
			if((int)lact % (int)((double)Constants.LACT / 5) == 0 || lact == 0 || lact == Constants.LACT)
			{
				if(lact > 3 || lact == 0) System.out.println((double)lact/(double)Constants.LACT*100);
				Utility.write_msg_file(22,(double)lact/(double)Constants.LACT*100);
			}
		}//lact
	
	
		System.out.println("\n BIO_PROBABILITY done ...");
		return bio;
	}
	/*** end BIO_PROBABILITY ***/
	
	


	/******************************************************************************
	CALVING_SEA: Finds sea when cow calved.
	- calving at start of month
	- e.g. -2 = October, 0 = December, 8 = August, 14 = February etc.
	******************************************************************************/
	//unsigned short calving_sea(unsigned short	sea,		//season at event
	//						   short			months)		//mim at event
	public int calving_sea(int mo)	
	{
		//	while (months > Constants.SEA) months -= Constants.SEA; 
		//	if(sea >= months) return sea - months + 1;
		//	else return sea - months + 13;
		if(mo <= 0) while (mo < 1) mo += Constants.SEA;
		if(mo >= 1) while (mo > Constants.SEA) mo -= Constants.SEA; 
		return mo;
	}
	/*** end CALVING_SEA ***/


	/******************************************************************************
	BIO_SMY: Calculates the level of milk production per curr.
	-  See BIO_PROBABILITY
	******************************************************************************/
	public BioRecord bio_smy(BioRecord bio)
	{
		int curr = 0;
		double sdev	= 0.0, xlm = 0.0, xum = 0.0, zxlm = 0.0, zxum = 0.0, pxlm = 0.0, pxum = 0.0;
		double yum[] = {0, 74,78,82,86,90,94,98,102,106,110,114,118,122,126,130};	//upper limits of levels
	
		System.out.println("\n BIO_SMY started ...");
		sdev = bio.getPa_cv();				//phenotypic coefficient of variation of daily of milk weights
		if(sdev == 0)						//test: no variability in milk production: all cows in class 8
		{
			bio.getSmy()[8] = 100;
			System.out.println("\n BIO_SMY done ...");
			return bio;
		}
		else if(sdev > 0.0)
		{
			//variability in milk production
			for(curr = 1; curr <= Constants.CURR; curr++)	//calculate predicted weight in every level
			{
				if(curr == 1) xlm = -9999;						//standardized lower limit of level m
				else xlm = (yum[curr-1] - 100)/sdev;
				xum  = (yum[curr] - 100)/sdev;					//standardized upper limit of level m
				if(curr == Constants.CURR) xum = 9999;
				zxlm = 1/Math.sqrt(2*Constants.PI)*Math.exp(-0.5*Math.pow(xlm,2));		//heigth of stand. normal distribution at point xlm
				zxum = 1/Math.sqrt(2*Constants.PI)*Math.exp(-0.5*Math.pow(xum,2));		//height of stand. normal distribution at point xum 
				pxlm = Utility.normal(xlm);								//cumulative probability left of xlm
				pxum = Utility.normal(xum);								//cumulative probability left of xum
				if((pxum - pxlm) > 0.000001) bio.getSmy()[curr] = 100 - sdev*(zxum - zxlm)/(pxum - pxlm);	//average production level (%)
				//	System.out.println("\ncurr:%2d xlm:%5.1f xum:%5.1f zxlm:%5.3f zxum:%5.3f pxlm:%7.3f pxum:%7.3f pxum-pxlm:%f smy:%7.3f",
				//	curr,xlm,xum,zxlm,zxum,pxlm,pxum,pxum - pxlm,bio.getSmy()[curr]);
				//	if(curr >= 2) System.out.println("  %f",bio.getSmy()[curr]-bio.getSmy()[curr-1]);
			}
		}
		//stop(0);
		System.out.println("\n BIO_SMY done ...");
		return bio;
	}
	/*** end BIO_SMY ***/

	
	
	
	/******************************************************************************
	PREDICT_305M: Determines the 305-d cummulative milk yield.
	- See notes in ESTIMATE_MILK()
	******************************************************************************/
	@SuppressWarnings("unused")
	public BioRecord predict_305m(BioRecord bio)
	{
		int curr = 0, lact = 0, mim = 0, sea = 0, msea = 0, i = 0;
	
		System.out.println("\n PREDICT_305M started ...");
		for(sea = 1; sea <= Constants.SEA; sea++)
			for(curr = 1; curr <= Constants.CURR; curr++)
				for(lact = 1; lact <= bio.getMaxlact(); lact++)
				{
				//	System.out.println("\n** sea:%d curr:%d lact:%d",sea,curr,lact);
					msea = sea;
					for(mim = 1; mim <= 10; mim++)	//10 months of 30.4 days = 305 days
					{
						bio.getM305()[curr][lact][sea] += bio.getMilkyield()[curr][lact][mim][0][msea];
				//		System.out.println(" %.1f",bio.getMilkyield()[curr][lact][mim][0][msea]/Constants.STAGE);
						msea = msea + 1; if(msea > Constants.SEA) msea = 1;
					}
				//	System.out.println(" :%.0f",bio.getm305[curr][lact][sea]);
				//	stop(0);
				}
		//print:
	//	for(lact = 1; lact <= Constants.LACT; lact++)
	//		for(curr = 1; curr <= Constants.CURR; curr++)
	//			for(sea = 1; sea <= Constants.SEA; sea++)
	//			{
	//				if(sea == 1) System.out.println("\n lact:%d curr:%2d ",lact,curr);
	//				System.out.println(" %5.0f",bio.getm305[curr][lact][sea]);
	//			}	
	//	stop(0);
		System.out.println("\n PREDICT_305M done ...");
		return bio;
	}
	/*** end PREDICT_305M ***/
	
	
	
	public TestRecord init_test_record(BioRecord bio)
	{
		int	lact = 0, curr = 0, mim = 0, mpreg = 0, sea = 0, w = 0;
		//test = (struct test_record*)malloc(sizeof(struct test_record));
		TestRecord test = new TestRecord();
		/*if(test == NULL)
		{
			System.out.println("\n INIT_TEST_RECORD: Insufficient memory available. ");
			stop(1);
		}
		else
		{*/
			for(w = 0; w <= 12; w = w + 12)
			for(sea = 1; sea <= Constants.SEA; sea++)
			for(curr = 1; curr <= Constants.CURR; curr++)
			for(lact = bio.getMaxlact(); lact >= 1; lact--)
				for(mim = Constants.MIM; mim >= 1; mim--)
					for(mpreg = 0; mpreg <= Constants.MPREG; mpreg++) if(mpreg == 0 || (mim - mpreg >= 2 && mim - mpreg <= bio.getLbm()))
					{
						test.getRpo()[curr][lact][mim][mpreg][sea+w] = -Constants.X;
						test.getPol()[curr][lact][mim][mpreg][sea+w] = Constants.POL;
					}
			System.out.println("\n INIT_TEST_RECORD done ...");
		/*}*/
		return test;
	}
	
	
	public double rpo_test(BioRecord bio, RpoRecord rpo, TestRecord test)
	{
		int	lact = 0, curr = 0, mim = 0, mpreg = 0, sea = 0, w = 12;
		double			rpodiff = 0.0, maxrpodiff = 0.0;
		
		for(sea = 1; sea <= Constants.SEA; sea++)
		for(curr = 1; curr <= Constants.CURR; curr++)
			for(lact = bio.getMaxlact(); lact >= 1; lact--)
				for(mim = Constants.MIM; mim >= 1; mim--)
					for(mpreg = 0; mpreg <= Constants.MPREG; mpreg++) if(mpreg == 0 || (mim - mpreg >= 2 && mim - mpreg <= bio.getLbm()))
					{
						if(rpo.getCow()[curr][lact][mim][mpreg][sea+w] != -Constants.X && test.getRpo()[curr][lact][mim][mpreg][sea+w] != -Constants.X)
						{
						//	System.out.println("\n %f %f",rpo.getCow()[curr][lact][mim][mpreg][sea+w],test.getRpo()[curr][lact][mim][mpreg][sea+w]);
							rpodiff = Math.abs(rpo.getCow()[curr][lact][mim][mpreg][sea+w] - test.getRpo()[curr][lact][mim][mpreg][sea+w]);
							if(rpodiff > maxrpodiff) maxrpodiff = rpodiff;
						}
						test.getRpo()[curr][lact][mim][mpreg][sea+w] = rpo.getCow()[curr][lact][mim][mpreg][sea+w];
					}
		return maxrpodiff;
	}
	
	
	
	@SuppressWarnings("unused")
	public int pol_test(BioRecord bio, PolRecord pol, TestRecord test)
	{
		int	lact = 0, curr = 0, mim = 0, mpreg = 0, sea = 0, w = 12;
		double			rpodiff = 0.0, maxrpodiff = 0.0;
		int				polcount = 0;
	
		for(sea = 1; sea <= Constants.SEA; sea++)
		for(curr = 1; curr <= Constants.CURR; curr++)
			for(lact = bio.getMaxlact(); lact >= 1; lact--)
				for(mim = Constants.MIM; mim >= 1; mim--)
					for(mpreg = 0; mpreg <= Constants.MPREG; mpreg++) if(mpreg == 0 || (mim - mpreg >= 2 && mim - mpreg <= bio.getLbm()))
					{
					//	if(curr == 8 && lact == 1)
					//		System.out.println("\n %d %d",pol.getCow()[curr][lact][mim][mpreg][sea+w],test.getpol[curr][lact][mim][mpreg][sea+w]);
						if(pol.getCow()[curr][lact][mim][mpreg][sea+w] != test.getPol()[curr][lact][mim][mpreg][sea+w])
							polcount++;
						test.getPol()[curr][lact][mim][mpreg][sea+w] = pol.getCow()[curr][lact][mim][mpreg][sea+w];
					}
		return polcount;
	}
	
		
	
	/******************************************************************************
	DP_POLICY: Calculates optimal insemination and replacement policy for individual
	cows, including the Retention Pay-Off (RPO = KEEP - REPLACE) for all possible states.
	- KEEP: involuntary replacement at end of stage.
	- REPLACE: animal is replaced at start of stage.
	- Repeat pure DP calculations untill results are stable. Then add the pol.getpolicy
	  to the front. 
	- dpp.getpolicy is the steady-state optimal policy.
	- pol.getpolicy is the user-defined policy. But pol.get[] == 8 is the part that is optimized by DP_POLICY
	  pol.getpolicy is only calculated after dpp.getpolicy has been determined.
	- Dynamic programming based on Van Arendonk[128]103 and Huirne[155]80. 
	  May 2003: milk class transition between months included (Houben[816]2986.
	- Repeatability milk production classes from Van Arendonk[128]119.
	- Net return due to genetic progress: 60% of 115 kg/year = $19.73 per year.
	  Add genetic progress to replacement value.
	- Assumes no temporary environmental effects (on milk production or reproduction).
	- Some erratic behavior may occur near conversion due to float rounding errors.
	
	//see Survival function for P(culling | cow survived until now) PUBH8420
	//RPO can never be more than (p_heifer - sell price cow) for average cow.
	
	DP_POLICY() and MARKOV_CHAIN() reconsiliation:   May 30, 2003
	-heifer0 should be the same as mkv.getCash()[0] when the same policies are followed.
	-debugging: start DP with just one stage and compare heifer0 with mkv.getCash()[0] with one
	            starting heifer and her replacements, for t=1 to 1
	       etc, add one stage to DP and compare heifer0 with mkv.getCash()[0] for t=1 to 2  etc.
				no more than 12 stages can be compared this way: if e.g. 15 stages were run, then
				DP would have 15 different stages, but MKV would have 12 different stages and
				3 stages doubled.  So the MKV policy would be different from the DP policy and
				therefore give different results.
				Make BIO_DEFAULT() as simple as possible.
	-check:     when DP has converged, then difference heifer0[t=1]-heifer0[t=2] equals 
	            mkv.getCash()[t] for one stage with steady state herd.
	******************************************************************************/
	@SuppressWarnings("unused")
	public int dp_policy(BioRecord bio, FutRecord fut, PolRecord pol, InsRecord ins, RpoRecord rpo, DelRecord del, MkvRecord mkv)
	{
		int 	futfile = 0;								//fut file: 0= do not use, 1 = use
		char	biofiletime[] = new char[100];
		char 	message[] = new char[100];
		
		int policy = bio.getPolicy();
		
		
		
		
		long			iter		= 0, counter = 0;		//iteration
		int				calc		= -1;		//calculation stage helper
		int				i = 0, j = 0, x = 0, k = 0, m = 0, t = 0, 
						w			= 12,		//0=year 1, 12=later years
						w1			= 12,
						lo			= 0,		//lowest level with meaningfull transition probabilities  
						hi			= 0,		//highes level with meaningfull transition probabilities
						next		= 0,		//phenotypic leven next lactation
						curr		= 0,		//phenotypic level current lactation
						lact		= 0,		//lactation number
						lact3		= 0,		
						mim			= 0,		//months in milk
						mpreg		= 0,		//months pregnant
						sea			= 0,		//season: Jan:1  Dec:12
						sea_		= 0,		//season help counter
						sea1		= 0,
						mpr			= 0,		//month of later pregnancy (Monsanto question)
						mop			= 0,		//month of later being open
						errors		= 0,		//calculation errors
						heif_itt	= 1,		//acceleration iteration
			//			help1[MAXU+1],			//pol array helper
			//			help2[MAXU+1],			//pol array helper
						heif_cnt	= 0,
						id			= 0,
						ii			= 0,		//swap helper
						jj			= 1;		//swap helper
		double			pa_progress = 0.0,		//weekly net return due to producing ability progress
				//		heifer0		= 0.0,		//heifer revenue stream at start of stage
						heifer1		= 0.0,		//heifer revenue stream at end of stage
						p_heifer	= 0.0,		//heifer price
						p_loss		= 0.0,		//$ loss due to involuntary culling
						p_fixed		= 0.0,		//fixed cost per stage
						ppreg		= 0.0,		//probability pregnancy
						reprocost	= 0.0,		//cost of heatdetection or fertility drugs (labor + drugs)
						breedcost	= 0.0,		//cost of inseminations (labor + semen/embryo)
						df			= 0.0,		//Utility.discount factor per week
						pheat		= 0.0,		//probability of heat detection
						p_sell0		= 0.0,		//cow salvage value start of stage (voluntary)
						p_sell1		= 0.0,		//cow salvage value end of stage (involuntary)
						p_sell2		= 0.0,		//cow salvage value end of stage (voluntary)
						keep		= 0.0,		//income stream if decision to keep (no serve)
						keep1		= 0.0,		//income stream of open cow (Monsanto question)
						serv		= 0.0,		//income stream if decision to serve
						serv1		= 0.0,
						repl		= 0.0,		//income stream if decision to replace
						buy			= 0.0,
						delay		= 0.0,
						pnext		= 0.0,
						yield		= 0.0,
						rev			= 0.0;
		double 			heif[] = new double[3+1];
		double			diff		= 0.0,
						pp			= 0.0,
						p1			= 0.0,
						p2			= 0.0,
						fut_in		= 0.0,
						fut_out     = 0.0,
						maxdiff		= 0.0, 
						sumfixed	= 0.0,		//sum of fixed cost -.get used to prevent too early convergence
						abort		= 0.0,
						pid			= 0.0;		//probability of involuntary culling
	
	//	FILE			*fp	= NULL;
	//	if((fp = fopen("dp.txt","a")) == NULL) { stop(1); }
	
		System.out.println("\n DP_POLICY started ... policy:"+bio.getPolicy());
		heif[0] = heif[1] = heif[2] = heif[3] = 0.0;
	//	if(bio.getpolicy >= 3) stop(1);
	//	System.out.println("\n policy: %d",bio.getpolicy);
	//	System.out.println("  startseason: %d",bio.getStart_sea());
	//	System.out.println("\n  serve: "); for(sea = 1; sea <= 2*Constants.SEA; sea++) System.out.println("%d ",bio.getserve[sea]);
	//	System.out.println("\n  enter: "); for(sea = 1; sea <= 2*Constants.SEA; sea++) System.out.println("%d ",bio.getenter[sea]);
	//stop(0);
		if(bio.getVwp()[1] < 2 || bio.getVwp()[2] < 2 
				|| bio.getVwp()[3] < 2 || bio.getLbm() > Constants.MIM-Constants.MPREG 
				|| bio.getVwp()[1] > bio.getLbm() || bio.getVwp()[2] > bio.getLbm() 
				|| bio.getVwp()[3] > bio.getLbm()) 
			{ System.out.println("\n bio.getvwp1:"+bio.getVwp()[1]+" or bio.getvwp2:"
				+bio.getVwp()[2]+" or bio.getvwp3:"+bio.getVwp()[3]+" or bio.getlbm:"
				+bio.getLbm()+" error"); Utility.stop(1); }
		
		/*System.out.println("bio.getVwp()[1] < 2: "+(bio.getVwp()[1] < 2));
		System.out.println("bio.getVwp()[2] < 2: "+(bio.getVwp()[2] < 2));
		System.out.println("bio.getVwp()[3] < 2: "+(bio.getVwp()[3] < 2));
		System.out.println("bio.getLbm() > Constants.MIM-Constants.MPREG: "+(bio.getLbm() > Constants.MIM-Constants.MPREG));
		
		System.out.println("bio.getVwp()[1] > bio.getLbm(): "+(bio.getVwp()[1] > bio.getLbm()));
		System.out.println("bio.getVwp()[2] > bio.getLbm(): "+(bio.getVwp()[2] > bio.getLbm()));
		System.out.println("bio.getVwp()[3] > bio.getLbm(): "+(bio.getVwp()[3] > bio.getLbm()));
		*/
		if(bio.getMaxmim() > Constants.MIM || bio.getMaxmim() <= bio.getLbm()) 
			{ System.out.println("\n bio.getmaxmim:"+bio.getMaxmim()+" or bio.getlbm:"
				+bio.getLbm()+" error"); Utility.stop(1); }
		
	//	reset_fut_record();												//resets fut record (if futfile == 0) then do not read fut file
		df = Utility.discount(bio.getInterest()/12,1);								//monthly Utility.discount factor
	//	bio.getpa_progress = 19.73 / 12;									//$ producing ability progress per month
		pa_progress = 0.0;
		p_loss = bio.getP_loss();											//$ loss from involuntary culling
		p_fixed = (bio.getP_fixed_labor() + bio.getP_fixed_other()) * Constants.STAGE;	//fixed cost for stage
		sea = bio.getStart_sea();
		do
		{
		//	hlp = fut.getCow()[0][1][1][0][ii];
		//	for(i = 1; i <= Constants.CURR; i++)
		//		for(j = 1; j <= Constants.LACT; j++)
		//			for(k = 1; k <= Constants.MIM; k++)
		//				for(m = 0; m <= Constants.MPREG; m++)
		//					fut.getCow()[i][j][k][m][ii] -= hlp;
		//System.out.println("\n fut0:%f",fut.getCow()[0][1][1][0][ii]); stop(0);
	
			counter++;											//do loop counter
			sea--; if(sea < 1) sea = Constants.SEA;						//season in previous stage for next iteration
			//calculation stage definition (replaced dpp_done and pol_done):
			//-2   : iterate until dpp steady state is reached, w = 12 (prices year 2 + later)
			//-1   : iterate until pol steady state is reached, w = 12 (prices year 2 + later)
			// 1-12: calculate remaining pol data, w = 12 (prices year 2 + later)
			//13-24: calculate remaining pol data, w = 0 (prices year 1)
			if(calc >= 0) calc++;
			if(calc >= 13) w = 0;
			if(sea == bio.getStart_sea()) iter++;					//year counter
			p_heifer = bio.getP_heifer()[sea+w];					//$ heifer price
			if(counter == 1) fut_in = fut.getEnter()[0]; else fut_in = fut_out;//future income difference before this stage
					//fut_in = fut.getenter[0] is ok when interest = 0%, but not when interest > 0%
	//*** ENTERING HEIFERS ******************************************************************
			pid = bio.getPind()[1][1][sea];
			delay = fut.getEnter()[ii] * df + pa_progress - p_fixed;//delay at start of stage until optimal time to enter
			heifer1 = fut.getEnter()[ii] * df + pa_progress;		//best decision one stage in future
			buy = - p_heifer + bio.getP_calf() - p_fixed;			//buy heifer at start of stage, keep at least 1 month
			if(sea == Constants.SEA) p_sell1 = bio.getSel()[1][1][1+w] * bio.getInvcullpct();
			else p_sell1 = bio.getSel()[1][1][sea+1+w] * bio.getInvcullpct();
			for(curr = 1; curr <= Constants.CURR; curr++)
			{
				buy += bio.getPnext()[0][0][0][0][curr] * (bio.getRev()[curr][1][1][0][sea+w]
					+ pid * (- p_loss + p_sell1 + heifer1) * df);//cull and replace at end of stage
				lo = bio.getLo()[curr][1][1][0];
				hi = bio.getHi()[curr][1][1][0];
				for(next = lo; next <= hi; next++)				//a heifer is always kept the first month	
				{
					buy += bio.getPnext()[0][0][0][0][curr] * bio.getPnext()[curr][1][1][0][next] 
						* (1 - pid) * fut.getCow()[next][1][2][0][ii] * df;
				}
			}
			if(bio.getEnter()[sea+w] == 1) delay = buy - 0.001;		//always buy heifer: delay < buy
			if(bio.getEnter()[sea+w] == 0) buy = delay - 0.001;		//always delay heifer: buy < delay
	
			if(counter == 13)									//forced buy: prevents early convergence with empty herd
				if(pol.getEnter()[13] == 0 && pol.getEnter()[14] == 0 && pol.getEnter()[15] == 0 && pol.getEnter()[16] == 0 && pol.getEnter()[17] == 0 && pol.getEnter()[18] == 0
				&& pol.getEnter()[19] == 0 && pol.getEnter()[20] == 0 && pol.getEnter()[21] == 0 && pol.getEnter()[22] == 0 && pol.getEnter()[23] == 0 && pol.getEnter()[24] == 0)
					delay = buy - 0.001;
	
			//fut.getHeifer0() = max2(buy,delay);
			fut.setHeifer0(Utility.max2(buy,delay));
			fut.getEnter()[jj] = fut.getHeifer0();
		//	System.out.println("\n\nbio.getenter:%d %d heifer1:%f  buy:%f  delay:%f  buy-delay:%f  heifer0:%f",bio.getenter[sea+w],sea,heifer1,buy,delay,buy-delay,fut.getheifer0);
		//	stop(0);
			if(buy < delay) pol.getEnter()[sea+w] = 0; else pol.getEnter()[sea+w] = 1;
		//	System.out.println("\n\n"); for(i = 1; i <= Constants.SEA; i++) System.out.println(" %d",bio.getenter[i]); for(i = 1; i <= Constants.SEA; i++) System.out.println(" %d",pol.getenter[i]);
		//	System.out.println("\n"); for(i = 1; i <= Constants.SEA; i++) System.out.println(" %d",pol.getenter[i]);
		//	System.out.println(" %d",pol.getenter[sea+w]);
			for(curr = 1; curr <= Constants.CURR; curr++)
			{
				for(lact = 1; lact <= bio.getMaxlact(); lact++)
				{
					if(lact <= 3) lact3 = lact; else lact3 = 3;
					for(mim = 1; mim <= Constants.MIM; mim++)
					{
	//*** OPEN COWS *************************************************************************
						mpreg = 0;
						pid = bio.getPind()[lact][mim][sea];
						p_sell0 = bio.getSel()[lact][mim][sea+w];					//salvage value start of stage
						if(sea == Constants.SEA) p_sell2 = bio.getSel()[lact][mim][1+w];		//salvage value end of stage
						if(sea <  Constants.SEA) p_sell2 = bio.getSel()[lact][mim][sea+1+w];
						p_sell1 = p_sell2 * bio.getInvcullpct();					//salvage value involuntary cull
						lo = bio.getLo()[curr][lact][mim][0];						//lowest level for transition probabilities
						hi = bio.getHi()[curr][lact][mim][0];						//highest level for transition probabilities
						if(mim >= bio.getMaxmim())									//last month open
						{	//an rpo is calculated for 1 to Constants.MIM: allows the calculation of the value of pregnancy
							//it is also a clean way of dealing with cows aborting after bio.getmaxmim and becoming open again
							serv = -Constants.X;
							keep = bio.getRev()[curr][lact][mim][0][sea+w] - p_fixed
								+ (1 - pid) * (p_sell2 + heifer1) * df			//end of stage (voluntary cull)
								+ pid * (- p_loss + p_sell1 + heifer1) * df;	//end of stage (involuntary cull)
							serv1 = serv;
							keep1 = keep;
						}
						else if(mim >= bio.getVwp()[lact3] && mim <= bio.getLbm())		//breeding opportunity
						{														//  at end of stage
							p1 = bio.getHeatdetectionrate()[mim][sea] * bio.getConceptionrate()[lact][mim][sea];
							p2 = bio.getHeatdetectionrate()[mim][sea];
							pheat = Utility.min2(p2 + (1 - p1) * p2 * 9.4 / 21,1.0);		//probility of service in a month
							ppreg = p1 + (1 - p1) * p1 * 9.4 / 21;				//21 days for pp, 30.4-21=9.4 days for (1-pp)*pp
							reprocost = bio.getReprocost()[mim][sea];				//cow/month
							breedcost = pheat * (bio.getP_semen() + bio.getTimebreeding() / 60 * bio.getP_labor());
	
						//old code while fixing the difference with markov_chain()
						//	p2 = bio.getheatdetectionrate[mim][sea] * (bio.getp_semen + bio.gettimebreeding / 60 * bio.getp_labor);
						//	breedcost = p2 + (1 - p1) * p2 * 9.4 / 21;			//1-p1 = eligible for estrus next time
						//	p1 = bio.getheatdetectionrate[mim][sea] * bio.getconceptionrate[lact][mim][sea];
						//	p2 = bio.getheatdetectionrate[mim][sea];
						//	pheat = min(p2 + (1 - p1) * p2 * 9.4 / 21,1.0);		
						//	ppreg = p1 + (1 - p1) * p1 * 9.4 / 21;				//probability of pregnancy in a month
						//	pelig = 1 + (1 - p1) * 9.4 / 21;					//eligible days in a month
	
							serv = bio.getRev()[curr][lact][mim][0][sea+w] - p_fixed - reprocost	
								 - (1 - pid) * breedcost * df					//try to breed, service after pid
								 + pid * (- p_loss + p_sell1 + heifer1) * df;
							serv1 = serv;
							for(next = lo; next <= hi; next++)
							{
								serv += bio.getPnext()[curr][lact][mim][0][next] * 
									((1 - pid) * (1 - ppreg) * fut.getCow()[next][lact][mim+1][0][ii]
									+ (1 - pid) * ppreg * fut.getCow()[next][lact][mim+1][1][ii]) * df;
							}
							keep = bio.getRev()[curr][lact][mim][0][sea+w] - p_fixed//do not try to breed
								+ pid * (- p_loss + p_sell1 + heifer1) * df;
							keep1 = keep;
							for(next = lo; next <= hi; next++)
							{
								keep += bio.getPnext()[curr][lact][mim][0][next]
									* (1 - pid) * fut.getCow()[next][lact][mim+1][0][ii] * df;
							}
							//constraint: always or never serv
							ins.getCow()[curr][lact][mim][sea+w] = serv - keep;		//insemination value at end of month
							//!!! IS INSEMINATION VALUE OVERESTIMATED? BECAUSE RISK OF PREGNANCY IS 30/21 Constants.X PREGRATE + TAKES ANOTHER MONTH FOR NEW INSEMINATION OPPORTUNITY
	
							if(bio.getServe()[sea+w] == 1) keep = serv - 0.001;			//always serv (RPO ok)
							if(bio.getServe()[sea+w] == 0) serv = keep - 0.001;			//never serv (RPO ok)
							//end: constraint
	
	//if(counter >= 12)
	//{
	//System.out.println("\nopen %d curr:%d lact:%d mim:%d mpreg:%d sea:%d",counter,curr,lact,mim,mpreg,sea);
	//System.out.println(" bio.getmilk:%.0f  cullmilk:%.0f",bio.getMilkyield()[curr][lact][mim][0][sea],bio.getcullmilk);
	//System.out.println("\n %f %f %f %f",bio.getRev()[curr][lact][mim][0][sea+w],reprocost,keep,serv);
	//System.out.println("\n %f %f %f %f",breedcost,pid,ppreg,heifer1);
	//System.out.println("\n %f %f %f %f",bio.getheatdetectionrate[mim][sea],bio.getp_semen,bio.getlaborcostbreeding,p_loss);
	//}
	
						}
						else if(mim <= bio.getVwp()[lact3] - 1 || mim >= bio.getLbm() + 1)	//no breeding opportunity
						{
							serv = -Constants.X;											//reset serv
							keep = bio.getRev()[curr][lact][mim][0][sea+w] - p_fixed
								+ pid * (- p_loss + p_sell1 + heifer1) * df;
							if(mim == 1) keep += bio.getP_calf();
							serv1 = serv;
							keep1 = keep;
	
	//if(mim == 1)
	//{
	//System.out.println("\nopen %d curr:%d lact:%d mim:%d mpreg:%d sea:%d",counter,curr,lact,mim,mpreg,sea);
	//System.out.println(" rev:%.0f  p_fixed:%.0f pid:%f p_loss:%.0f p_sell1:%.0f heifer1:%.0f",
	//	   bio.getRev()[curr][lact][mim][0][sea+w],p_fixed,pid,p_loss,p_sell1,heifer1);
	//}
	
	
							for(next = lo; next <= hi; next++)
							{
								keep += bio.getPnext()[curr][lact][mim][0][next]
									* (1 - pid) * fut.getCow()[next][lact][mim+1][0][ii] * df;
	//if(mim == 1)
	//System.out.println("\n next:%d  keep:%.0f",next,keep);
	
							}
	
	if(keep <= serv)
	{
		System.out.println("\nERROR open "+counter+" curr:"+curr+" lact:"+lact+" mim:"+mim+" mpreg:"+mpreg+" sea:"+sea+" serv:"+serv+" >= keep:"+keep);
		Utility.stop(0);
	}
	
	
	
						}
						else { System.out.println("\n OPENerror"); Utility.stop(1);}
						repl = p_sell0 + fut.getHeifer0();
						if(mim == 1) repl += bio.getP_calf();
						///////////////////////////////////////////////////////////////////////////
						//Monsanto question:
						//what is the economic value of getting a cow pregnant today vs 'mpr' months later?
						//mpr = month in milk at start of pregnancy (begin of month)
						if(mim >= bio.getVwp()[lact3] && mim <= bio.getLbm()) for(mpr = mim+1; mpr <= bio.getLbm()+1; mpr++)
						{	//cow gets pregnant "mpr" months after mim, begin of month
							pp = 0.0; for(next = lo; next <= hi; next++)
							{
								pp += bio.getPnext()[curr][lact][mim][0][next]
									* (1 - pid) * fut.getPreg()[next][lact][mim+1][mpr][ii] * df;
							}
							fut.getPreg()[curr][lact][mim][mpr][jj] = Utility.max2(serv1,keep1)+pp;
	
	
	//if(curr == 8 && lact == 1 && calc == 24)
	//{
	//System.out.println("\npreg %d curr:%d lact:%d mim:%d",counter,curr,lact,mim);
	//System.out.println(" mpr:%d  fut=max(serv1:%.2f keep1:%.2f)+pp:%.2f = %.2f",mpr,serv1,keep1,pp,fut.getPreg()[curr][lact][mim][mpr][jj]);
	//if(mim == bio.getlbm) stop(0);
	//}
	
	
						}
						///////////////////////////////////////////////////////////////////////////
						//reverse of Monsanto question:
						//what is the economic loss of optimal decision now vs. staying open and optimal decision "mop" months later
						//mop = month in milk at start of pregnancy (begin of month)
	
						//	if(bio.getserve[sea+w] == 1) keep1 = serv1 - 1;			//always serv (RPO ok)
						//	if(bio.getserve[sea+w] == 0) serv1 = keep1 - 1;			//never serv (RPO ok)
	
						if(mim < bio.getMaxmim()) for(mop = mim+1; mop <= bio.getMaxmim(); mop++)		
						{	//cow remains open "mop" months after mim, begin of month
							pp = 0.0; for(next = lo; next <= hi; next++)
							{
								pp += bio.getPnext()[curr][lact][mim][0][next]
								* (1 - pid) * fut.getOpen()[next][lact][mim+1][mop][ii] * df;
							}
							fut.getOpen()[curr][lact][mim][mop][jj] = Utility.max2(serv1,keep1)+pp;
	
	//if(curr == 8 && lact == 1 && calc == 24)
	//{
	//System.out.println("\nopen %d curr:%d lact:%d mim:%d ",counter,curr,lact,mim);
	//System.out.println(" mop:%d  fut=max(serv1:%.2f keep1:%.2f) + pp:%.2f = %.2f",mop,serv1,keep1,pp,fut.getOpen()[curr][lact][mim][mop][jj]);
	//if(mim == bio.getlbm) stop(0);
	//}
	
	
						}
						///////////////////////////////////////////////////////////////////////////////
	
	
						yield = bio.getMilkyield()[curr][lact][mim][0][Constants.SEA] / Constants.STAGE;
						rev = bio.getRev()[curr][lact][mim][0][Constants.SEA] / Constants.STAGE;
						if((yield < bio.getCullmilk() && mim >= 4)								//too low milk production
						|| (rev < bio.getCullmargin() && mim >= 4))								//too low margin
						{																	//policy: replace
							fut.getCow()[curr][lact][mim][0][jj] = repl;
							pol.getCow()[curr][lact][mim][0][sea+w] = 0;
						//	rpo.getCow()[curr][lact][mim][0][sea+w] = -XXXX;						//initialized to -Constants.X: rpo does not exist
						//	System.out.println("\nopen:%d sea:%d curr:%d lact:%d mim:%d",counter,sea,curr,lact,mim);
						//	System.out.println(" yield:%.1f <? cullyield:%.1f, or",yield,bio.getcullmilk);
						//	System.out.println(" rev:%.1f <? cullmargin:%.1f",rev,bio.getcullmargin);
						//	System.out.println("\a");
						//	stop(0);
						}
						else if(policy == 1 || policy == 2)									//policy: OPTIMAL or KEEPPREG
						{
							fut.getCow()[curr][lact][mim][0][jj] = Utility.max3(serv,keep,repl);
							pol.getCow()[curr][lact][mim][0][sea+w] = Utility.max3u(repl,keep,serv) - 1;
						//	rpo.getCow()[curr][lact][mim][0][sea+w] = max2(serv,keep) - repl;
	
	if(mim <= bio.getVwp()[lact3] - 1 || mim >= bio.getLbm() + 1) if(pol.getCow()[curr][lact][mim][0][sea+w] == 2)
	{
		System.out.println("\nERROR open "+counter+" curr:"+curr+" lact:"+lact+" mim:"+mim+" mpreg:"+mpreg+" sea:"+sea+" serv:"+serv+" keep:"+keep);
		System.out.println(" pol:"+pol.getCow()[curr][lact][mim][0][sea+w]);
	//	stop(0);
	}
	
	
	//{
	//System.out.println("\n\n error: open %d curr:%d lact:%d mim:%d mpreg:%d sea:%d  pol:%d",counter,curr,lact,mim,mpreg,sea,
	//	   pol.getCow()[curr][lact][mim][mpreg][sea+w]);
	//System.out.println("\n policy:%d repl:%f keep:%f serv:%f  pol:%d",policy,repl,keep,serv,max3u(repl,keep,serv) - 1);
	//}
	
	
						}
						else if(policy == 3)												//policy: KEEPALL
						{
							fut.getCow()[curr][lact][mim][0][jj] = Utility.max2(serv,keep);
							pol.getCow()[curr][lact][mim][0][sea+w] = Utility.max2u(keep,serv);
						//	rpo.getCow()[curr][lact][mim][0][sea+w] = max2(serv,keep) - repl;
						}
						rpo.getCow()[curr][lact][mim][0][sea+w] = Utility.max2(serv,keep) - repl;
						fut.getOpen()[curr][lact][mim][mim][jj] = fut.getCow()[curr][lact][mim][0][jj];
	
	
	//if((pol.getCow()[curr][lact][mim][mpreg][sea+w] >= 1 && rpo.getCow()[curr][lact][mim][mpreg][sea+w] < 0)
	//|| (pol.getCow()[curr][lact][mim][mpreg][sea+w] == 0 && rpo.getCow()[curr][lact][mim][mpreg][sea+w] > 0))
	//{
	//System.out.println("\nerror(pol rpo): open %d curr:%d lact:%d mim:%d mpreg:%d sea:%d  pol:%d rpo:%.0f mlk:%.2f",counter,curr,lact,mim,mpreg,sea,
	//	   pol.getCow()[curr][lact][mim][mpreg][sea+w],rpo.getCow()[curr][lact][mim][mpreg][sea+w],bio.getMilkyield()[curr][lact][mim][mpreg][sea]/Constants.STAGE);
	//System.out.println("\npolicy:%d repl:%f keep:%f serv:%f  pol:%d",policy,repl,keep,serv,max3u(repl,keep,serv) - 1);
	//stop(0);
	//}
	
	//*** PREGNANT COWS *********************************************************************
						for(mpreg = 1; mpreg <= Constants.MPREG; mpreg++) if(mim - mpreg >= bio.getVwp()[lact3] && mim - mpreg <= bio.getLbm())
						{
							p_sell0 = bio.getSel()[lact][mim][sea+w];
							if(sea == Constants.SEA) p_sell2 = bio.getSel()[lact][mim][1+w];//salvage value end of stage
							if(sea <  Constants.SEA) p_sell2 = bio.getSel()[lact][mim][sea+1+w];
							p_sell1 = p_sell2 * bio.getInvcullpct();
							abort = bio.getAbort()[mpreg];						//risk of abortion at end of month 1 to 8
							if(mpreg <= Constants.MPREG - 1)							//pregnant, not calving
							{
								lo = bio.getLo()[curr][lact][mim][0];			//lowest level for transition probabilities
								hi = bio.getHi()[curr][lact][mim][0];			//highest level for transition probabilities
								keep = bio.getRev()[curr][lact][mim][mpreg][sea+w] - p_fixed
									 + pid * (- p_loss + p_sell1 + heifer1) * df;
								for(next = lo; next <= hi; next++)
								{
									//System.out.println("curr: "+curr+" lact: "+lact+" mim: "+mim+" next: "+next);
									double tempPNext = bio.getPnext()[curr][lact][mim][0][next];
									
									//System.out.println("next: "+next+" lact: "+lact+" mim+1: "+(mim+1)+" mpreg+1: "+(mpreg+1)+" ii: "+ii);
									double tempCow1 = fut.getCow()[next][lact][mim+1][mpreg+1][ii];
									double tempCow2 = fut.getCow()[next][lact][mim+1][0][ii];
									
									keep += tempPNext * (1 - pid)
											  * ((1 - abort) * tempCow1
											  + abort * tempCow2) * df;
									
									/*keep += bio.getPnext()[curr][lact][mim][0][next] * (1 - pid)
										  * ((1 - abort) * fut.getCow()[next][lact][mim+1][mpreg+1][ii]
										  + abort * fut.getCow()[next][lact][mim+1][0][ii]) * df;*/
								}
							}
							else if(mpreg == Constants.MPREG && lact <= bio.getMaxlact() - 1)	//pregnant, calving
							{
								lo = bio.getLo()[curr][lact][mim][sea];			//transition to next lactation
								hi = bio.getHi()[curr][lact][mim][sea];			//transition to next lactation
								keep = bio.getRev()[curr][lact][mim][mpreg][sea+w] - p_fixed
									 + pid * (- p_loss + p_sell1 + heifer1) * df;
	
	//help = 0.0;
								for(next = lo; next <= hi; next++)
								{	
									keep += bio.getPnext()[curr][lact][mim][sea][next]
										 * (1 - pid) * fut.getCow()[next][lact+1][1][0][ii] * df;
	
	//help += bio.getpnext[curr][lact][0][next] * (1 - pid) * fut.getCow()[next][lact+1][1][0][ii] * df;
	
								}
	
	//if(keep > 0) System.out.println("\n curr:%d lact:%d mim:%d  help:%f keep:%f",curr,lact,mim,help,keep);
	
	
							//	stop(0);
							}
							else if(mpreg == Constants.MPREG && lact == bio.getMaxlact())			//end of last lactation
							{
	//old
	//							keep = bio.getRev()[curr][lact][mim][mpreg][sea+w] - p_fixed
	//								 + pid * (- p_loss + p_sell1 + heifer1) * df	//involuntary cull
	//								 + (1 - pid) * (p_sell2 + heifer1) * df;		//voluntary cull
	//old
	
	
	//new: bio.getMaxlact()
	//  -.get all future lactations are the same, so no forced culling in last lactation
	//  -.get markov chain cash/year in steady state may be different from dp cash/year
								lo = bio.getLo()[curr][lact][mim][sea];	
								hi = bio.getHi()[curr][lact][mim][sea];		
								keep = bio.getRev()[curr][bio.getMaxlact()][mim][mpreg][sea+w] - p_fixed
									 + pid * (- p_loss + p_sell1 + heifer1) * df;
								for(next = lo; next <= hi; next++)
								{	
									keep += bio.getPnext()[curr][lact][mim][sea][next]
										 * (1 - pid) * fut.getCow()[next][lact][1][0][ii] * df;
								}
	//new
							}
							else { System.out.println("\nPREG error"); Utility.stop(1); }
							repl = p_sell0 + fut.getHeifer0();
							if(policy == 1)											//policy: OPTIMAL
							{
								fut.getCow()[curr][lact][mim][mpreg][jj] = Utility.max2(keep,repl);
								pol.getCow()[curr][lact][mim][mpreg][sea+w] = Utility.max2u(repl,keep) - 1;
							}
							else if(policy == 2 || policy == 3)						//policy: KEEPPREG or KEEPALL
							{
								fut.getCow()[curr][lact][mim][mpreg][jj] = keep;
								pol.getCow()[curr][lact][mim][mpreg][sea+w] = 1;
							}
							if(mpreg == 1) fut.getPreg()[curr][lact][mim][mim][jj] = fut.getCow()[curr][lact][mim][1][jj];		  
							rpo.getCow()[curr][lact][mim][mpreg][sea+w] = keep - repl;
	
	
	
	//if((pol.getCow()[curr][lact][mim][mpreg][sea+w] >= 1 && rpo.getCow()[curr][lact][mim][mpreg][sea+w] < 0)
	//|| (pol.getCow()[curr][lact][mim][mpreg][sea+w] == 0 && rpo.getCow()[curr][lact][mim][mpreg][sea+w] > 0))
	//{
	//System.out.println("\nerror(pol rpo): preg %d curr:%d lact:%d mim:%d mpreg:%d sea:%d  pol:%d rpo:%.0f mlk:%.2f",counter,curr,lact,mim,mpreg,sea,
	//	   pol.getCow()[curr][lact][mim][mpreg][sea+w],rpo.getCow()[curr][lact][mim][mpreg][sea+w],bio.getMilkyield()[curr][lact][mim][mpreg][sea]/Constants.STAGE);
	//stop(0);
	//}
	
	
	//System.out.println("\npreg %d curr:%d lact:%d mim:%d mpreg:%d sea:%d  pol:%d  rpo:%.0f  mlk:%.2f",counter,curr,lact,mim,mpreg,sea,
	//	   pol.getCow()[curr][lact][mim][mpreg][sea+w],rpo.getCow()[curr][lact][mim][0][sea+w],bio.getMilkyield()[curr][lact][mim][mpreg][sea]/Constants.STAGE);
	
						}//end for(mpreg) PREGNANT COWS
					}//end for(mim)
				}//end for(lact)
				if(counter == 1 && (curr == 1 || curr == 4 || curr == 7 || curr == 10 || curr == 15))
				{
					System.out.println((double)curr/Constants.CURR*100);
					Utility.write_msg_file(4,(double)curr/Constants.CURR*100);
				}
			}//end for(curr)
	//*** DP BOTTOM *************************************************************************
			fut_out = fut.getHeifer0();					//future revenue begin of stage
			for(i = Constants.SEA*2; i >= 2; i--) fut.getCash()[i] = fut.getCash()[i-1];
			fut.getCash()[1] = fut_out - df * fut_in;
	
			// accelerate convergence ////////////////////////////////////////////////////////
			// method: fit exponential distribution; CHESS-RO, Ruud Huirne, Wageningen University
			if(sea == 1 && counter >= 36 && df < 1.0 && calc == -1 && heif_itt <= 6) heif[++heif_cnt] = fut.getHeifer0();
			if(sea == 1 && heif_cnt == 3             && calc == -1 && heif_itt <= 6) //set heif_itt for the # accelerations
			{
			//	System.out.println("\n [%d] %.2f %.2f %.2f",heif_cnt,heif[1],heif[2],heif[3]);
				System.out.println(" A");
				heif[0] = estimate_heifer0(heif[1],heif[2],heif[3]);
				heif[1] = heif[2] = heif[3] = 0.0;
				heif_cnt = 0;
				heif_itt++;
				if(heif[0] != -Constants.X && Math.abs(heif[0] - fut_out) > 20)
				{
					System.out.println(heif[0]);
					diff = heif[0] - fut_out;
					fut_out += diff;
					fut.getEnter()[jj] += diff;
					//fut.getHeifer0() += diff;
					fut.setHeifer0(fut.getHeifer0()+diff);
					for(curr = 1; curr <= Constants.CURR; curr++)
						for(lact = 1; lact <= Constants.LACT; lact++)
							for(mim = 1; mim <= Constants.MIM; mim++)
								for(mpreg = 0; mpreg <= Constants.MPREG; mpreg++) if(mpreg == 0 || (mim - mpreg >= 2 && mim - mpreg <= 15))
									fut.getCow()[curr][lact][mim][mpreg][jj] += diff;	
				}
			}
			//note at convergence when df > 0: 
			//   df = 1 - (fut.getCash()[1] / fut_out);
			//   fut.getCash()[1] = fut_out - df * fut_out 
			//   fut.getCash()[1] = (1 - df) * fut_out;
			//   fut.getCash()[1] = fut_out - fut_in;
			//--------------------------------------.get this is the same solution as policy iteration!
			/////////////////////////////////////////////////////////////////////////////////////////
	
			//June/July 2005:
			//Equivalent annual annuity = EAA = NPV / pvifa, see Keown et al. Foundations of finance, page 554
			//Calculate EAA for last 12 months and compare to previous 12 months, if same: converge
			//See 'Equivalent period costs.xls' in Meadows review
			//fut_cash[1] can change, but the change is perhaps due to a different path the algorithm takes and 
			//the real cash change is far in the future.  So say fut_cash is always 100 but then jumps to 101: this 
			//may be a result of a different path with 100 cash now and 101 cash somewhere (far) in the future. Calculating
			//EAA for the last 12 stages does not do justice to that and is not necessarily correct.
			//EAA calculated for the same 12 stages in MARKOV_CHAIN_STATS is correct and may be different from the EAA calculated
			//here in DP_POLICY.  Under stead-state conditions, the EAA in both the DP_POLICY and MARKOV_CHAIN_STATS should
			//be the same.
			mkv.setCash12(0.0); for(i = 1; i <= Constants.SEA; i++) mkv.setCash12(mkv.getCash12() + fut.getCash()[i]);
			//Because EAA is DP_POLICY is not correct, this codes does not matter anymore:
			//pvifa = 0.0; if(counter < Constants.SEA) for(i = 1; i <= counter; i++) pvifa = pvifa + pow(df,i-1);	//present-value interest factor for an annuity
			//else for(i = 1; i <= Constants.SEA; i++) pvifa = pvifa + pow(df,i-1);	//present-value interest factor for an annuity
			//System.out.println(" df=%f pvifa=%f",df,pvifa);
	
			maxdiff = 0.0;
			for(i = 1; i <= Constants.SEA; i++) 
			{
				maxdiff += Math.abs(fut.getCash()[i] - fut.getCash()[i+12]);
			//	System.out.println("\nfut_cash[%d]:%.2f  fut_cash[%d]:%.2f  maxdiff%.2f",i,fut_cash[i],i+12,fut_cash[i+12],maxdiff);
			}
			sumfixed -= p_fixed;
			if(fut_out == sumfixed) maxdiff = 1.0;//prevents early convergence when alwaysbuy == 0 
		//	System.out.println("\n%d calc:%2d w:%2d sea:%2d fut_out:%.0f fut_cash:%6.2f  yr_cash:%.3f  opp:%6.2f  maxdiff:%f ",
		//		counter,calc,w,sea,fut_out,fut_cash[1],cash12,buy-delay,maxdiff);	
		//	for(i = 1; i <= 12; i++) System.out.println("%d",pol.getenter[i+w]);
	
			System.out.println("\n "+counter+" calc:"+calc+" w:"+w+" sea:"+sea+" fut:"+fut_out+" cash:"+fut.getCash()[1]+" yr:"+mkv.getCash12()+" conv:"+maxdiff+" pol:"+pol.getEnter()[sea+w]);
		//	System.out.println(" %5.2f %5.2f %5.2f %5.2f",fut_in,fut.getCash()[0],fut.getCash()[1],fut.getCash()[2]);
	
	//		if(calc > 0)
		//	fSystem.out.println(fp,"\n %d calc: %d w: %d sea: %d df: %f fut_out: %f cash: %f yr: %f conv: %f heifer0: %f heifer1: %f fut_in: %f fut_cash1: %f fut_cash2: %f  heif[0]: %f",
		//	counter,calc,w,sea,df,fut_out,fut.getCash()[1],cash12,maxdiff,fut.getheifer0,heifer1,fut_in,fut.getCash()[1],fut.getCash()[2],heif[0]);
	
	
		//	if(sea == bio.getStart_sea()) System.out.println(".");
			if(sea == bio.getStart_sea() && counter >= 2) Utility.write_msg_file(4,maxdiff);
			//////////////////////////////////////////////////////////////////////////////////////////////////////
			//###########################
	/*
			if(calc >= 1)
			{
				if(calc == 1) read_cow_file();
				for(id = 1; id <= maxid; id++) if((double)cow[id].getDim() / 30.4 >= calc)
				{
					mim = (unsigned short)ceil((double)cow[id].getDim() / 30.4);
	
	System.out.println("\n id:%d dim:%d mim:%d",id,cow[id].getDim(),mim);
	
					dp_cow(policy,calc,sea,ii,id,mim,fut.getheifer0,heifer1,p_loss,p_calf,p_fixed,df);
	
	stop(1);
				}
	
	
	//			if(calc == 24) doe nog een keer met left + find inbetween estimates
				if(calc == 24) write_cow_file();
			}
	*/
			//###########################################
			////////////////////////////////////////////////////////////////////////////////////////////////////////
	
	
	
			if(calc >= 1)													//last 2 years, store pregnancy value
			{
				mkv.getOpportunitycost()[sea+w] = buy - delay;
				for(curr = 1; curr <= Constants.CURR; curr++)
				{
					for(lact = 1; lact <= bio.getMaxlact(); lact++)
					{
						if (lact <= 3) lact3 = lact; else lact3 = 3;
						for(mim = 1; mim <= Constants.MIM; mim++)
						{
							for(mpreg = 0; mpreg <= Constants.MPREG; mpreg++)
							{
								if(mpreg == 0 && mim >= bio.getVwp()[lact3]+1 && mim <= bio.getLbm()+1)	//open cows
								{	//fut.getCow()[jj] is at start of stage
									rpo.getPrg()[curr][lact][mim][0][sea+w] =
										fut.getCow()[curr][lact][mim][1][jj] - fut.getCow()[curr][lact][mim][0][jj];
								//	System.out.println("\n curr:%2d lact:%2d mim:%2d sea+w:%2d",curr,lact,mim,sea+w);
								//	System.out.println(" fut:%5.0f - fut:%5.0f = prg:%4.0f",fut.getCow()[curr][lact][mim][1][jj],fut.getCow()[curr][lact][mim][0][jj],rpo.getprg[curr][lact][mim][0][sea+w]);
								//	System.out.println(" rpo[0]:%4.0f",rpo.getCow()[curr][lact][mim][0][sea+w]);
								//	System.out.println(" rpo[1]:%4.0f",rpo.getCow()[curr][lact][mim][1][sea+w]);
								//	System.out.println(" ins:%3.0f",ins.getCow()[curr][lact][mim][sea+w]);
								//	System.out.println(" pol:%d",pol.getCow()[curr][lact][mim][0][sea+w]);
	
							//		//can insemination value have different sign than pregnancy value?  YES
							//		//especially late lactation cows have ins < 0 and prg > 0, probably because prg
							//		//is after the cost of breeding, and ins is before the cost of breeding (with semen cost etc).
							//		if((ins.getCow()[curr][lact][mim-1][sea+w] > 0 && rpo.getprg[curr][lact][mim][0][sea+w] < 0)
							//		|| (ins.getCow()[curr][lact][mim-1][sea+w] < 0 && rpo.getprg[curr][lact][mim][0][sea+w] > 0))
							//		{	//-.getif ins < 0 then sometimes prg > 0
							//			System.out.println("\nSIGN! curr:%2d lact:%2d mim:%2d sea+w:%2d",curr,lact,mim,sea+w);
							//			System.out.println(" ins:%8.3f  prg:%8.3f  pol:%d",ins.getCow()[curr][lact][mim-1][sea+w],rpo.getprg[curr][lact][mim][0][sea+w],
							//				pol.getCow()[curr][lact][mim-1][0][sea+w]);
							//			if(ins.getCow()[curr][lact][mim-1][sea+w] < 0 && rpo.getprg[curr][lact][mim][0][sea+w] > 0) System.out.println(" *");
							//			if(ins.getCow()[curr][lact][mim-1][sea+w] > 0 && rpo.getprg[curr][lact][mim][0][sea+w] < 0) stop(0);
							//		}
							//		if(rpo.getprg[curr][lact][mim][0][sea+w] < -0.1)
							//		{	//-.getif prg < 0 then always ins < 0
							//			System.out.println("\nPRG<0 curr:%2d lact:%2d mim:%2d sea+w:%2d",curr,lact,mim,sea+w);
							//			System.out.println(" ins:%8.3f  prg:%8.3f  pol:%d",ins.getCow()[curr][lact][mim-1][sea+w],rpo.getprg[curr][lact][mim][0][sea+w],
							//				pol.getCow()[curr][lact][mim-1][0][sea+w]);
							//			if(ins.getCow()[curr][lact][mim-1][sea+w] > 0) stop(1);
							//		}
							//	
	
								}
	
								if(mpreg >= 1 && mim - mpreg >= bio.getVwp()[lact3] && mim - mpreg <= bio.getLbm())	//pregnant cows
								{
									rpo.getPrg()[curr][lact][mim][mpreg][sea+w] =
										fut.getCow()[curr][lact][mim][mpreg][jj] - fut.getCow()[curr][lact][mim][0][jj];
								//	System.out.println("\n mim:%d mpreg:%d sea:%d",mim,mpreg,sea+w);
								//	System.out.println(" fut:%.0f - fut:%.0f = prg:%.0f",fut.getCow()[curr][lact][mim][mpreg][jj],fut.getCow()[curr][lact][mim][0][jj],rpo.getprg[curr][lact][mim][mpreg][sea+w]);
								//	if(rpo.getprg[curr][lact][mim][mpreg][sea+w] < -1000) stop(0);
								}
							}
							//cost of longer open: economics of optimal decision for open cows now vs not getting pregnant
							//   but keeping open cows "mop" months longer
							for(mop = mim; mop <= Constants.MIM; mop++)
							{
								del.getOpen()[curr][lact][mim][mop][sea+w] = 
									 fut.getOpen()[curr][lact][mim][mop][jj] - fut.getCow()[curr][lact][mim][0][jj];
								if(curr == 8 && lact == 1 && calc == 24)
								{
							//		System.out.println("\nopen curr:%2d lact:%2d mim:%2d mop:%d sea+w:%2d",curr,lact,mim,mop,sea+w);
								//	System.out.println(" pol:%d",pol.getCow()[curr][lact][mim][0][sea+w]);
							//		System.out.println(" ins:%3.0f",ins.getCow()[curr][lact][mim][sea+w]);
								//	System.out.println(" rpo(open):%5.0f",rpo.getCow()[curr][lact][mim][0][sea+w]);
								//	System.out.println(" rpo(preg):%5.0f",rpo.getCow()[curr][lact][mim][1][sea+w]);
								//	System.out.println("  %.2f",rpo.getCow()[curr][lact][mim][1][sea+w]-rpo.getCow()[curr][lact][mim][0][sea+w]);
							//		System.out.println(" del.getOpen():%4.0f",del.getOpen()[curr][lact][mim][mop][sea+w]);
							//		System.out.println(" fut:%5.0f - fut:%5.0f = del.getOpen():%4.0f",fut.getOpen()[curr][lact][mim][mop][jj],fut.getCow()[curr][lact][mim][0][jj],del.getOpen()[curr][lact][mim][mop][sea+w]);
							//		if(mop == Constants.MIM) stop(0);
								}
							}
							//Monsanto question: economics of getting open cow pregnant now vs "mpr" months later
							if(mim >= bio.getVwp()[lact3]+1 && mim <= bio.getLbm()) for(mpr = mim+1; mpr <= bio.getLbm()+1; mpr++) 
							{
								del.getPreg()[curr][lact][mim][mpr][sea+w] = 
									 fut.getPreg()[curr][lact][mim][mpr][jj] - fut.getCow()[curr][lact][mim][1][jj];
							//	if(curr == 8 && lact == 1 && calc == 12)
							//	{
							//		System.out.println("\npreg curr:%2d lact:%2d mim:%2d mpr:%d sea+w:%2d",curr,lact,mim,mpr,sea+w);
							//		System.out.println(" pol:%d",pol.getCow()[curr][lact][mim][0][sea+w]);
							//		System.out.println(" rpo:%5.0f",rpo.getCow()[curr][lact][mim][1][sea+w]);
							//		System.out.println(" fut.getPreg():%5.0f - fut.getCow():%5.0f = del.getPreg():%4.0f",fut.getPreg()[curr][lact][mim][mpr][jj],fut.getCow()[curr][lact][mim][1][jj],del.getPreg()[curr][lact][mim][mpr][sea+w]);
							//		if(mpr == bio.getlbm+1) stop(0);
							//	}
							}
	
						}//end mim
					}//end lact
				}//end curr
			}//end calc >= 1
	
	
			if(calc == -1 && sea == bio.getStart_sea() && maxdiff < Constants.MAXDIFF && counter >= 48)	//pol converged
			{
				System.out.println("*");
				calc = 0;													//now go run calc= 1-24 (w=12 and w=0)
				for(i = 1; i <= Constants.SEA*2; i++) pol.getEnter()[i] = bio.getEnter()[i];	//reset pol.getenter[] record, will be filled again
			}
			if(calc == 12 && futfile == 1) Utility.write_warm_file(jj);				//save warm.txt for warm start
			if(calc == 24)
			{
	
	// needs more checking and reconcile with cull rate:   1/cullrate x 12 = rkeep of a heifer entering the herd, in steady state
				//remaining breeding period	(rbreed) = how much longer can this open cow optimally be bred (idea: David Galligan)
				//remaining open period (rlact) = how much longer can this open cow optimally be kept in the current lactation = how long is RPO(open) > 0?
				//remaining herd life (rpl) = how much longer is this cow expected to remain in the herd (my idea)
	
				for(lact = bio.getMaxlact(); lact >= 1; lact--)
				{
					if(lact <= 3) lact3 = lact; else lact3 = 3;
					for(mim = Constants.MIM; mim >= 1; mim--)
					{
						for(mpreg = 0; mpreg <= Constants.MPREG; mpreg++) if(mpreg == 0 || (mim - mpreg >= bio.getVwp()[lact3] && mim - mpreg <= bio.getLbm()))	
						{
							for(i = bio.getStart_sea(); i <= bio.getStart_sea() + 23; i++)
							{
								sea = i; 
								if(sea >= bio.getStart_sea() + Constants.SEA) while(sea > Constants.SEA*2) sea = sea - Constants.SEA;
								else while(sea > Constants.SEA) sea = sea - Constants.SEA;
								sea1 = i + 1;
								if(sea1 >= bio.getStart_sea() + Constants.SEA) while(sea1 > Constants.SEA*2) sea1 = sea1 - Constants.SEA;
								else while(sea1 > Constants.SEA) sea1 = sea1 - Constants.SEA;
								sea_ = sea; while(sea_ > Constants.SEA) sea_ = sea_ - Constants.SEA;
								pid = bio.getPind()[lact][mim][sea_];
						//		pid = 0.0;
	
								for(curr = 1; curr <= Constants.CURR; curr++)
								{
									if(pol.getCow()[curr][lact][mim][mpreg][sea] >= 1) rpo.getRpl()[curr][lact][mim][mpreg][sea] = 1.0;
	
	//if(i == 1 && curr == 8) if(mpreg == 0) System.out.println("\n");
	//if(i == 1 && curr == 8) System.out.println("\ni:%2d curr:%2d lact:%2d mim:%2d mpreg:%2d sea_:%2d sea:%2d sea1:%2d pid:%.3f",i,curr,lact,mim,mpreg,sea_,sea,sea1,pid);
	//if(i == 1 && curr == 8) System.out.println(" pol:%d ",pol.getCow()[curr][lact][mim][mpreg][sea]);
	
	pp = 0.0;
	
									if(pol.getCow()[curr][lact][mim][mpreg][sea] >= 1) 
									{
										if(mpreg == 0 && mim < Constants.MIM)								//not pregnant
										{
											lo = bio.getLo()[curr][lact][mim][0];					//lowest level for transition probabilities
											hi = bio.getHi()[curr][lact][mim][0];					//highest level for transition probabilities
											if(pol.getCow()[curr][lact][mim][0][sea] <= 1)			//cow is not inseminated
												for(next = lo; next <= hi; next++)
												{
													rpo.getRpl()[curr][lact][mim][0][sea] += 
														bio.getPnext()[curr][lact][mim][0][next] * (1 - pid) * rpo.getRpl()[next][lact][mim+1][0][sea1];
	
	
	//if(mim == 0 && curr == 8) System.out.println("\n %f",del.getrkeep[curr][lact][mim][0][sea]);
	
													pp += bio.getPnext()[curr][lact][mim][0][next];
	
												}
	
											else if(pol.getCow()[curr][lact][mim][0][sea] == 2)		//cow is inseminated
											{
												p1 = bio.getHeatdetectionrate()[mim][sea_] * bio.getConceptionrate()[lact][mim][sea_];
												ppreg = p1 + (1 - p1) * p1 * 9.4 / 21;			//21 days for pp, 30.4-21=9.4 days for (1-pp)*pp
		
	//ppreg = 1.0;
	
												for(next = lo; next <= hi; next++)
												{
													rpo.getRpl()[curr][lact][mim][0][sea] += 
														bio.getPnext()[curr][lact][mim][0][next] * (1 - pid) *
															((1 - ppreg) * rpo.getRpl()[next][lact][mim+1][0][sea1] + ppreg * rpo.getRpl()[next][lact][mim+1][1][sea1]);
	
												pp += bio.getPnext()[curr][lact][mim][0][next];
	
												}
											}
	//if(i == 1 && curr == 8) System.out.println("*");
	
										}
										else if(mpreg >= 1 && mpreg <= Constants.MPREG - 1)				//pregnant, not calving
										{
											abort = bio.getAbort()[mpreg];	
											lo = bio.getLo()[curr][lact][mim][0];					//lowest level for transition probabilities
											hi = bio.getHi()[curr][lact][mim][0];					//highest level for transition probabilities
											for(next = lo; next <= hi; next++)
											{
												rpo.getRpl()[curr][lact][mim][mpreg][sea] += 
													bio.getPnext()[curr][lact][mim][0][next] * (1 - pid) * 
													((1 - abort) * rpo.getRpl()[next][lact][mim+1][mpreg+1][sea1] + 
													 (abort) * rpo.getRpl()[next][lact][mim+1][0][sea1]);
	
											pp += bio.getPnext()[curr][lact][mim][0][next];
	
											}
	//if(i == 1 && curr == 8) System.out.println("$");
	
										}
										else if(mpreg == Constants.MPREG && lact <= bio.getMaxlact() - 1)				//pregnant, calving
										{
											lo = bio.getLo()[curr][lact][mim][sea_];				//lowest level for transition probabilities
											hi = bio.getHi()[curr][lact][mim][sea_];				//highest level for transition probabilities
											for(next = lo; next <= hi; next++)
											{
												rpo.getRpl()[curr][lact][mim][mpreg][sea] += 
													bio.getPnext()[curr][lact][mim][sea_][next] * (1 - pid) * rpo.getRpl()[next][lact+1][1][0][sea1];	
	
											pp += bio.getPnext()[curr][lact][mim][sea_][next];
											}
	//if(i == 1 && curr == 8) System.out.println("#");
	
	
										}
										else
										{
	//if(i == 1 && curr == 8) System.out.println("&");
	
										}
									}
		
	//if(i == 1 && curr == 8) System.out.println(" lo:%2d  hi:%2d  %.3f rkeep:%5.2f",lo,hi,pp,rpo.getrpl[curr][lact][mim][mpreg][sea]);
	//if(i == 1 && curr == 8) if (mpreg == 9) stop(0);
	//if(i == bio.getStart_sea() && pol.getCow()[curr][lact][mim][mpreg][sea] >= 1) stop(0);
	
								}
							}
						}
					}
				}
	//stop(0);
	
	
	/*old, but seemed ok
				//rkeep > rbreed when delayed breeding is optimal.
				for(i = bio.getStart_sea(); i <= bio.getStart_sea() + Constants.SEA * 2; i++)
					for(lact = Constants.LACT; lact >= 1; lact--)
						for(mim = Constants.MIM; mim >= 1; mim--)
							for(mpreg = 0; mpreg <= Constants.MPREG; mpreg++)
								for(curr = 1; curr <= Constants.CURR; curr++)
								{
									sea_ = i + mim - 1;
									if(sea_ > bio.getStart_sea() + Constants.SEA) w = 12; else w = 0;
									sea = sea_; while(sea > Constants.SEA) sea = sea - Constants.SEA;
									sea1 = sea + 1; if(sea1 > Constants.SEA) sea1 = 1;
									if(sea_ >= bio.getStart_sea() + Constants.SEA) w1 = 12; else w1 = 0;	
	
									pid = bio.getPind()[lact][mim][sea+w];
	
	System.out.println("\nstart:%d i:%d curr:%2d lact:%d mim:%d sea_:%d w:%d sea:%d pid:%.3f",bio.getStart_sea(),i,curr,lact,mim,sea_,w,sea,pid);
	System.out.println(" pol:%d ",pol.getCow()[curr][lact][mim][0][sea+w]);
	
									if(pol.getCow()[curr][lact][mim][0][sea+w] >= 1) del.getropen[curr][lact][mim][sea+w] = 1.0;
									if(pol.getCow()[curr][lact][mim][0][sea+w] == 2) del.getrbreed[curr][lact][mim][sea+w] = 1.0;
									lo = bio.getlo[curr][lact][mim][0];						//lowest level for transition probabilities
									hi = bio.gethi[curr][lact][mim][0];						//highest level for transition probabilities
									for(next = lo; next <= hi; next++)
									{
										del.getropen[curr][lact][mim][sea+w] += 
											bio.getpnext[curr][lact][mim][0][next] * (1 - pid) * del.getropen[next][lact][mim+1][sea1+w1];
										del.getrbreed[curr][lact][mim][sea+w] += 
											bio.getpnext[curr][lact][mim][0][next] * (1 - pid) * del.getrbreed[next][lact][mim+1][sea1+w1];
									}
	
	System.out.println(" lo:%2d  hi:%2d  rbreed:%5.2f  ropen:%5.2f",lo,hi,del.getrbreed[curr][lact][mim][sea+w],del.getropen[curr][lact][mim][sea+w]);
	if(mim >= bio.getvwp && mim <= bio.getlbm) System.out.println("*");
	if (curr == 15) stop(0);
								}
	*/
	
			}//calc == 24
			jj = ii; if(ii == 0) ii = 1; else ii = 0;	//swap ii and jj
			if(counter == 10000 || (counter == 200 && fut_out == sumfixed))
			{
				Utility.strcpy(message,"DP_POLICY not converged, biofiletime = ");
				Utility.strcat(message,biofiletime);
				//System.out.println("\a");
				Utility.system(message);
				Utility.write_warm_file(jj);	
				Utility.write_msg_file(5,0);
				return 0;	//not converged
			}
	
		} while(calc != 24);
		System.out.println("\n DP_POLICY done ... "+(double)(counter-36)/12+" years  calc: "+calc);
		for(i = 1; i <= 24; i++) System.out.println(pol.getEnter()[i+w]);
		
		Utility.write_msg_file(6,mkv.getCash12());
	//	stop(1);
	//	fclose(fp);
		return 1;			//converged
	}
	/*** end DP_POLICY ***/
	
	
	
	
	/******************************************************************************
	ESTIMATE_HEIFER0: Predict fut.getheifer0 when it converges when df > 0.
	- Same as procedure calc_enrgpd in unit calcdp1.pas (Ruud Huirne, CHESS-RO, Wageningen University)
	- Fits a negative exponential function to the heifer0 estimates.
	******************************************************************************/
	public double estimate_heifer0(double a1, double a2, double a3)
	{
		double	param_elab = 0.0, param_a = 0.0, param_b = 0.0, a1a2 = 0.0, a2a3 = 0.0;
	
		a1a2 = a1 - a2;
		a2a3 = a2 - a3;
	//	System.out.println("\n a1:%.2f a2:%f a3:%f a1a2:%f a2a3:%f",a1,a2,a3,a1a2,a2a3);
		if(a2a3 != 0 && (a1a2 * a2a3 - a2a3 * a2a3) != 0) 
		{
			param_elab = a1a2 / a2a3;
			param_b = a1a2 * a1a2 * a1a2 / (a1a2 * a2a3 - a2a3 * a2a3);
		//	System.out.println("\n param_elab:%f  param_b:%f ",param_elab,param_b);
			if(param_elab > 1.1)	//if param_elab > 1 then the change slows over time and convergence is possible
			{
				param_a = a1 - param_b / param_elab;
			//	System.out.println(" param_a:%f ",param_a);
			//	stop(0);
				return param_a;
			}
			else
			{
				System.out.println(":ratio: "+param_elab);
				return -Constants.X;
			}
		}
		else return -Constants.X;
	}
	/*** end ESTIMATE_HEIFER0 ***/
	

	
	
	/******************************************************************************
	ESTIMATE_MILK: Predict the future milk production in the lactation.
	- Does not include effect of pregnancy.
	- Estimate 305-d milk production, given transition probabilities within lactation.
	- See also PREDICT_305M()
	******************************************************************************/
	@SuppressWarnings("unused")
	void estimate_milk(BioRecord bio)
	{
		int	curr = 0, lact = 0, mim = 0, sea = 0, i = 0, j = 0, m = 0, next = 0,
						lo = 0, hi = 0, msea = 0;
		double			mlk = 0.0, sum = 0.0;
	
		System.out.println("\n ESTIMATE_MILK started ...");
		EstRecord est = new EstRecord(); 
		est = InitFunctions.init_est_record(est);							//initialized and resets est record
		for(curr = 1; curr <= Constants.CURR; curr++)
			for(lact = 1; lact <= 3; lact++)
			{
			//	System.out.println("\n curr:%d lact:%d",curr,lact);
				for(mim = 1; mim <= Constants.MIM; mim++)
				{
					est.getProb()[curr][lact][mim][mim][curr] = 1.0;
					for(m = mim + 1; m <= Constants.MIM; m++)
					{
						sum = 0.0;
						for(i = 1; i <= Constants.CURR; i++)
							for(j = 1; j <= Constants.CURR; j++)
							{
								est.getProb()[curr][lact][mim][m][j] += bio.getPnext()[i][lact][m-1][0][j] * est.getProb()[curr][lact][mim][m-1][i];
	//if(lact == 1 && curr == 7) 
	//System.out.println("\n curr:%d lact:%d mim:%d m:%d i:%d j:%d  pnext:%.6f  prob:%.6f",curr,lact,mim,m,i,j,bio.getpnext[i][lact][m-1][j],est.getprob[curr][lact][mim][m][j]);
	//if(lact == 2) 
	//stop(0);						
							}
	//if(lact == 1 && curr == 7) 
	//stop(0);
					}
				}
			}
	
		for(sea = 1; sea <= Constants.SEA; sea++)
			for(lact = 1; lact <= 3; lact++)
			{
			//	System.out.println("\n* sea:%d lact:%d",sea,lact);
				for(mim = 1; mim <= Constants.MIM; mim++)
				{
					msea = sea;						
					for(m = mim; m <= Constants.MIM; m++)		//sea applies to mim, not to m!
					{
						for(curr = 1; curr <= Constants.CURR; curr++)
						{
							for(i = 1; i <= Constants.CURR; i++)
							{
								est.getMlk()[curr][lact][mim][sea][m] += est.getProb()[curr][lact][mim][m][i] * bio.getMilkyield()[i][lact][m][0][msea];
						//		System.out.println("\n curr:%d lact:%d mim:%d sea:%d m:%d msea:%d i:%d mlk:%.0f  %.0f",curr,lact,mim,sea,m,msea,i,est.getmlk[curr][lact][mim][sea][m],bio.getMilkyield()[i][lact][m][0][msea]/Constants.STAGE);
							}
						}
						msea = msea + 1; if(msea > Constants.SEA) msea = 1;
					}
				}
	//stop(0);
			}
	
	
	/* Predict 305-d milk below seems correct. I use 305-d milk to adjust the transition probabilities between lactations
	   to adjust for days open in previous lactation (Van Arendonk[674]169), but that happens in bio_probability() 
	   while bio_probability() is also needed to calculate the 305-d milk here. So bio_probability would have to be called 2x.
	   This is too complex: instead first use predict_305m() to predict 305-day milk without transition probabilities witin
	   lactation, then use bio_probability()
		//predict 305-d milk yield:
		for(sea = 1; sea <= Constants.SEA; sea++)
			for(curr = 1; curr <= Constants.CURR; curr++)
				for(lact = 1; lact <= 3; lact++)
				{
				//	System.out.println("\n** sea:%d curr:%d lact:%d",sea,curr,lact);
					msea = sea;
					for(m = 1; m <= 10; m++)	//10 months of 30.4 days = 305 days
					{
						bio.getm305[curr][lact][sea] += est.getmlk[curr][lact][1][sea][m];
					//	System.out.println(" %.1f",est.getmlk[curr][lact][1][sea][m]/Constants.STAGE);
						if(lact == 3) for(i = 4; i <= Constants.LACT; i++) bio.getm305[curr][i][sea] = bio.getm305[curr][3][sea];
						msea = msea + 1; if(msea > Constants.SEA) msea = 1;
					}
				//	System.out.println(" :%.0f",bio.getm305[curr][lact][sea]);
				//	stop(0);
				}
		for(lact = 1; lact <= 3; lact++)
			for(curr = 6; curr <= 10; curr++)
				for(sea = 1; sea <= Constants.SEA; sea++)
				{
					if(sea == 1) System.out.println("\n lact:%d curr:%2d ",lact,curr);
					System.out.println(" %5.0f",bio.getm305[curr][lact][sea]);
				}	
	//	stop(0);
	*/
	
	
	/*
		curr = 6;
		lact = 1;
		mim = 1;
		sea = 1;
		for(i = 1; i <= Constants.CURR; i++)
		{
			System.out.println("\n");
			for(m = 1; m <= Constants.MIM; m++)		
				System.out.println(" %3.0f",est.getprob[curr][lact][mim][m][i]*100);
		}
		System.out.println("\n\n MILK curr:%d",curr);
		for(m = 1; m <= Constants.MIM; m++)		
			System.out.println(" %.0f",est.getmlk[curr][lact][mim][sea][m]/Constants.STAGE);
		System.out.println("  305:%.0f",bio.getm305[curr][lact][sea]);
	
		System.out.println("\n MILK curr:8");
		for(m = 1; m <= Constants.MIM; m++)		
			System.out.println(" %.0f",est.getmlk[8][lact][mim][sea][m]/Constants.STAGE);
		System.out.println("  305:%.0f",bio.getm305[8][lact][sea]);
	
		stop(0);
	*/
	
		est = InitFunctions.free_est_record(est);
		System.out.println("\n ESTIMATE_MILK done ...");		
		return;
	}/*** end ESTIMATE_MILK ***/
	
	
	
	
	/******************************************************************************
	MARKOV_CHAIN: Calculates the fraction cows and results in a certain state based
	on the policy.  Calculates the probability of entering a state in the next stage.
	- Entering heifers are calculated at the end of every stage.
	- For how = 0: simulate starting herd/cow over time for LP, no youngstock enters
	- For how = 1: keepfrac + replfrac == 1.0
	- For how = 2, keepfrac == 1.0 if excessive heifers are entered.
	
	- For how = 1, the stage is started with < 1.0 cows (non heifers) due to involuntary
	  culling at the end of the previous stage.  Then mkv.getcow7110 are added at the 
	  end of the stage to make sure a fraction of 1.0 cows started the stage and were
	  either kept or replace. At the end of the stage the involuntary culling takes
	  place. Therefore, at the end of the stage the fraction cows present < 1.0.
	- For how = 2, the stage is also started with < 1.0 cows (non heifers) due to involuntary
	  culling at the end of the previous stage. Then the fraction entering heifers is
	  added to calculate how many voluntary replacements are needed to start the stage
	  with a fraction 1.0 (after voluntary replacements are culled at start of stage).
	  At the end of the stage the involuntary culling takes place. Therefore, at the
	  end of the stage the fraction cows present < 1.0.
	******************************************************************************/
	@SuppressWarnings("unused")
	public void markov_chain(int pool,			//how is mkv.getheiferpool determined (entering heifers):
					  int maxstages,	//number of stages to simulate
					  int all,			//calculate all results
					  BioRecord bio,
					  MkvRecord mkv,
					  PolRecord pol,
					  RpoRecord rpo,
					  InsRecord ins,
					  DelRecord del)
	{	//pool:0 = LP simulate starting herd/cow over time, no replacements
		//pool:1 = DP unlimited heifers, sell calves, ignore starting youngstock
		//pool:2 = GA determines #entering + #leaving (based on ranked RPO)
		//pool:9 = simulate starting herd over time, heiferpool detemines #entering, no replacements
		int	t = 0, lact = 0, curr = 0, mim = 0, mpreg = 0, sea = 0, next = 0, cs0 = 0, cs1 = 0, dry = 0,
						decision = 0, ii = 1, jj = 0, repeat = 0, age = 0, lo = 0, hi = 0, w = 0, lact3 = 0, sea1 = 0,
						alldelayed = 0;
		int	i,j,k,m;
		double			mkvsum = 0.0, cost = 0.0, pheat1 = 0.0;
		double			pstate = 0.0, pid = 0.0, ppreg = 0.0, p_sell0 = 0.0, p_sell1 = 0.0, p_loss = 0.0, p_sell2 = 0.0,
						df0 = 1.0, df1 = 1.0, result = 0.0, reprocost = 0.0, breedcost = 0.0,
						p_heifer = 0.0, abort = 0.0, pfs = 0.0, pp = 0.0, pheat = 0.0, delig = 0.0, hlp = 0.0,
						pvw = 0.0, p1 = 0.0, p2 = 0.0;
	
	//	FILE			*fp	= NULL; 
	//	if((fp = fopen("mkv.txt","a")) == NULL) { stop(1); }
	
	//	System.out.println("\n MARKOV_CHAIN started ... startseason:%d  maxstages:%d all:%d",bio.getStart_sea(),maxstages,all);
		System.out.println("\n MARKOV_CHAIN ...");
		if(all == 0) mkv = InitFunctions.reset_mkv_record(mkv,0,maxstages);
		if(all == 0) mkv = InitFunctions.reset_mkv_record(mkv,1,maxstages);
		if(all == 1) mkv = InitFunctions.reset_mkv_record(mkv,3,maxstages);	//reset mkv-record, except mkv.getCow()[][][][][2] and mkv.getFs()[][][][2] and mkv.getHeifer()[][][2]
	
		//reset_mkv_record(3);
	
	//	mkvsum = mkv.getheiferpool[1] = 1.0;
	//	for(t = 2; t <= MAXU; t++) mkv.getheiferpool[t] = 0.0266;
	//	for(t = 1; t <= MAXU; t++) mkv.getheiferpool[t] = 0.05;
	//	System.out.println(" heiferpool from CHROM_TO_HEIFERPOOL()");
	//	for(t = 1; t <= 20; t++) System.out.println("\n mkv.getheiferpool[%d]:%f",t, mkv.getheiferpool[t]);
	//  stop(0);
	
		sea = bio.getStart_sea();					//season at t = 1
		for(i = 0; i <= Constants.CURR; i++)				//load starting cows + heifers in [ii]
			for(j = 0; j <= bio.getMaxlact(); j++)	//  starting herd is stored in [2]
				for(k = 0; k <= Constants.MIM; k++)
					for(m = 0; m <= Constants.MPREG; m++)
					{
						mkv.getCow()[i][j][k][m][ii] = mkv.getCow()[i][j][k][m][2];		//cows in states
						if(m == 0) mkv.getFs()[i][j][k][ii] = mkv.getFs()[i][j][k][2];	//cows in states eligible for first service
						if(m == 0) mkv.getVw()[i][j][k][ii] = mkv.getVw()[i][j][k][2];	//cows in states eligible for vwp
					}
		for(i = 0; i <= Constants.AGE; i++)
			for(m = 0; m <= Constants.MPREG; m++)
				mkv.getHeifer()[i][m][ii] = mkv.getHeifer()[i][m][2];					//heifers in states
	
	
		//if all delayed -.get no cows in steady state
		if(bio.getAnalysis() == 1)	//steady-state
		{
			alldelayed = 1;
			for(t = 1; t <= maxstages; t++)
				if(pol.getEnter()[t] == 1) alldelayed = 0;
		}
	
		p_loss = bio.getP_loss();
		if(all == 1) Utility.timetaken(3);
		for(t = 1; t <= maxstages; t++)
		{
			if(t >= 13) w = 12; else w = 0;
			sea1 = sea + 1; if(sea + 1 == 13) sea1 = 1;
	
			p_heifer = bio.getP_heifer()[0] + bio.getP_heifer()[sea+w];					//heifer price
		//	df0 = Utility.discount(bio.getInterest()/12,(double)t-1);						//discounting done in MARKOV_CHAIN_STATS()
		//	df1 = Utility.discount(bio.getInterest()/12,(double)t);							//discounting done in MARKOV_CHAIN_STATS()
	
		//	if(t % 20 == 0) System.out.println(" %d",t);
		//	mkvsum = 0.0;
		//	for(i = 0; i <= Constants.CURR; i++)
		//		for(j = 0; j <= Constants.LACT; j++)
		//			for(k = 0; k <= Constants.MIM; k++)
		//				for(m = 0; m <= Constants.MPREG; m++)
		//					mkvsum += mkv.getCow()[i][j][k][m][ii];	//entering stage
		//	System.out.println("\nt:%d mkvsum:%f",t,mkvsum);
	
		//	mkvsum = 0.0;
		//	for(i = 0; i <= Constants.AGE; i++)
		//		for(m = 0; m <= Constants.MPREG; m++)
		//			mkvsum += mkv.getHeifer()[i][m][ii];
		//	System.out.println(" mkvsum:%f",mkvsum);
		//	stop(0);
	
	/****** MOVE HEIFERS THROUGH THE Constants.STAGE ***************************************/
			for(age = 0; age <= Constants.AGE-1; age++)					//move female calves through time:		
				for(mpreg = 0; mpreg <= Constants.MPREG; mpreg++)
				{
					pstate = mkv.getHeifer()[age][mpreg][ii];
					if(all == 1) mkv.getFcalfcost()[t] += pstate * p_heifer / 23;//age 2 - 24 months, 1st month bottom
				//	if(mpreg == Constants.MPREG)							//pregnant, calving this month
				//		mkv.getheiferpool[t] += pstate;			//  see bottom:heifers entering
					if(mpreg >= 1 && mpreg <= Constants.MPREG-1)			//pregnant, not calving this month
						mkv.getHeifer()[age+1][mpreg+1][jj] += pstate;
					else if(age >= 15 && mpreg == 0)			//not pregnant, gets pregnant now
						mkv.getHeifer()[age+1][1][jj] += pstate;
					else if(mpreg == 0)							//not pregnant, stays open 1 more month
						mkv.getHeifer()[age+1][0][jj] += pstate;
				}
		//	for(age = 0; age <= Constants.AGE; age++)
		//		for(mpreg = 0; mpreg <= Constants.MPREG; mpreg++)
		//		{
		//			if(age == 0 && mpreg == 0) System.out.println("\n");
		//			if(mpreg == 0) System.out.println("%2d ",age);
		//			System.out.println("%.4f ",mkv.getHeifer()[age][mpreg][ii]);
		//			if(mpreg == Constants.MPREG) System.out.println("\n");
		//		}
		//	stop(0);
	
	
	/****** MOVE COWS THROUGH THE Constants.STAGE ******************************************/
			for(lact = 1; lact <= bio.getMaxlact(); lact++)				//move cows through time:	
			{
				if(t == 1 && all == 1 && (lact == 1 || lact == 3 || lact == 5 || lact == 7  || lact == 9 || lact == bio.getMaxlact())) 
				{
					System.out.println((double)lact/bio.getMaxlact()*50);	
					Utility.write_msg_file(15,(double)lact/bio.getMaxlact()*50);
				}
	
				if(lact <= 3) lact3 = lact; else lact3 = 3;
				for(curr = 1; curr <= Constants.CURR; curr++)
				{
					for(mim = 1; mim <= Constants.MIM; mim++)
					{
	//*** OPEN COWS *************************************************************************
						mpreg = 0;
						pstate = mkv.getCow()[curr][lact][mim][0][ii];
						if(mim <= bio.getVwp()[lact3]) pfs = pstate; else pfs = mkv.getFs()[curr][lact][mim][ii];
						if(mim <= bio.getVwp()[lact3]) pvw = pstate; else pvw = mkv.getVw()[curr][lact][mim][ii];
						if(pfs > pstate) { System.out.println("\nMARKOV_CHAIN error t:"+t+" curr:"+curr+" lact:"+lact+" mim:"+mim+" pstate:"+pstate+" < pfs:"+pfs); Utility.stop(1); }
						if(pstate > 0.00000000000001 && alldelayed == 0)	//if pstate == 0 then no need to calculate
						{
							decision = pol.getCow()[curr][lact][mim][0][sea+w];
							if(decision == Constants.POL) { System.out.println("\n decision error open t="+t); Utility.stop(1); }
							if(decision != 0 && decision != 1 && decision != 2) { System.out.println("\n decision error open = "+decision); Utility.stop(1); }
	
	//if(decision == 2)
	//{
	//	System.out.println("\nMKV t:%d  curr:%d  lact:%d  mim:%d  pstate:%f %d %d  sea:%d",
	//		t,curr,lact,mim,pstate,pol.getCow()[curr][lact][mim][0][t],dpp.getCow()[curr][lact][mim][0][sea],sea);
	//	stop(0);
	//}
	
							pid = bio.getPind()[lact][mim][sea];
							if(bio.getMilkyield()[curr][lact][mim][0][sea] >= 0.01) dry = 0; else dry = 1;
							p_sell0 = bio.getSel()[lact][mim][sea+w];
							if(sea == Constants.SEA) p_sell2 = bio.getSel()[lact][mim][1+w];//salvage value end of stage
							if(sea <  Constants.SEA) p_sell2 = bio.getSel()[lact][mim][sea+1+w];
							p_sell1 = p_sell2 * bio.getInvcullpct();
							//-p_heifer + p_calf1 had to be added in this stage to make
							//  fitness == fut_out in DP_POLICY
							lo = bio.getLo()[curr][lact][mim][0];					//lowest meaningful milk class transition
							hi = bio.getHi()[curr][lact][mim][0];					//highest meaningful milk class transition
							if(all == 1)
							{
								if(rpo.getCow()[curr][lact][mim][0][sea+w] != -Constants.X) mkv.getRpo()[t] += pstate * rpo.getCow()[curr][lact][mim][0][sea+w];
								if(rpo.getCow()[curr][lact][mim][0][sea+w] != -Constants.X) mkv.getRpo_cnt()[t] += pstate;
								if(rpo.getCow()[curr][lact][mim][0][sea+w] > 0) mkv.getValue()[t] += pstate * (rpo.getCow()[curr][lact][mim][0][sea+w] + p_sell0);
								else mkv.getValue()[t] += pstate * p_sell0;
								mkv.getValue_cnt()[t] += pstate;
	
	
	
	
	//if(mkv.getValue_cnt()[t] > 0)
	//System.out.println("\n curr:%d lact:%d mim:%d mpreg:%d sea:%d rpo:%.0f p_sell0:%.0f  value:%.0f  cnt:%f  %.0f",
	//	   curr,lact,mim,mpreg,sea,rpo.getCow()[curr][lact][mim][0][sea+w],p_sell0,mkv.getValue()[t],mkv.getValue_cnt()[t],mkv.getValue()[t] / mkv.getValue_cnt()[t]);
	//if(mim == Constants.MIM) stop(0);
	
							}
							if(mim >= bio.getMaxmim())										
							{
								if(decision == 1)								//decision: keep
								{
									mkv.getKeep()[curr][t] += pstate;					//finishes stage
									if(all == 1)
									{
										mkv.getOpen()[curr][t] += pstate;
										mkv.getPind()[curr][t] += pstate * pid;			//all cows culled
										mkv.getRepl()[curr][t] += pstate * (1 - pid);	//all cows culled
										mkv.getCash()[t] += pstate * (bio.getRev()[curr][lact][mim][0][sea+w] * df0
												 	 + pid * (- p_loss + p_sell1) * df1
													 + (1 - pid) * p_sell2 * df1);	//all cows are replaced at end of stage
										mkv.getRev()[t] += pstate * bio.getRev()[curr][lact][mim][0][sea+w] * df0;
										mkv.getMilkyield()[mim][t] += pstate * bio.getMilkyield()[curr][lact][mim][0][sea];
										mkv.getMilkyield_cnt()[mim][t] += pstate;
										if(dry == 1) mkv.getFeedcost()[0][t] += pstate * bio.getDmi()[curr][lact][mim][0][sea] * bio.getP_tmrdmdry() * df0;
										if(dry == 0) mkv.getFeedcost()[1][t] += pstate * bio.getDmi()[curr][lact][mim][0][sea] * bio.getP_tmrdmlac() * df0;
										if(dry == 0) mkv.getDmi()[mim][t] += pstate * bio.getDmi()[curr][lact][mim][0][sea];
										if(dry == 0) mkv.getDmi_cnt()[mim][t] += pstate;
										if(dry == 0) mkv.getMilking()[t] += pstate;
										mkv.getCowsales()[t] += pstate * pid * p_sell1 * df1;		//involuntary cull, end of state
										mkv.getCowsales()[t] += pstate * (1 - pid) * p_sell2 * df1;	//voluntary cull, end of state
										mkv.getLoss()[t] += pstate * pid * p_loss * df1;//loss due to involuntary culling
										mkv.getLact()[lact][t] += pstate;
										mkv.getDim()[t] += pstate * ((double)mim * Constants.STAGE - (Constants.STAGE / 2));
										mkv.getLaborcost()[t] += pstate * bio.getTimecows()[curr][lact][mim][0][sea] / 60 * bio.getP_labor() * df0;
										mkv.getLabortime()[t] += pstate * bio.getTimecows()[curr][lact][mim][0][sea];
										mkv.getVetcost()[t] += pstate * bio.getVetcost()[mim] * df0;
										mkv.getPind_mim()[mim][t] += pstate * pid;
										mkv.getRepl_mim()[mim][t] += pstate * (1 - pid);
										mkv.getStatus()[mim][0] += pstate;
										mkv.getDysopen()[t] += pstate * (double)(mim);
										mkv.getOthercost()[t] += pstate * bio.getP_other_lact() * Constants.STAGE * df0;
										mkv.getHz_open()[mim][0] += pstate;				//enter state
									//	mkv.getHz_open()[mim][1] += pstate * (1 - pid);	//repl open: should be at start of month (small error)
										mkv.getHz_open()[mim][2] += pstate * pid;		//pind open: at end of month
										//this is not quite right, as only pid is involuntary culled, but it is difficult to calculate in markov_chain_stats()
										//furthermore, it only affects the hazard in mim > bio.getmaxmim which is in a way "in"voluntary culling
										if(rpo.getCow()[curr][lact][mim][0][sea+w] != -Constants.X) mkv.getPind_rpo()[curr][lact][mim][0] += pstate * pid * rpo.getCow()[curr][lact][mim][0][sea+w];
										if(rpo.getCow()[curr][lact][mim][0][sea+w] != -Constants.X) mkv.getPind_rpo_cnt()[curr][lact][mim][0] += pstate * pid;
										//do we ever keep open cows that do not pay for themselves?
										if(bio.getRev()[curr][lact][mim][0][sea+w] < 0.0) mkv.getKeep_rev0()[t] += pstate;
									}
								}
								else if(decision == 0)							//decision: replace
								{
									if(all == 1)
									{
										mkv.getRepl()[curr][t] += pstate;				//does not finish stage
										mkv.getRepl_mim()[mim][t] += pstate;
										mkv.getCash()[t] += pstate * p_sell0 * df0;
										mkv.getCowsales()[t] += pstate * p_sell0 * df0;	//all cows sold at begin of stage
										mkv.getHz_open()[mim][0] += pstate;				//enter state
										mkv.getHz_open()[mim][1] += pstate;				//repl open: should be at start of month (small error)
	
									//	mkv.getRepl_open_rev()[curr][mim][t] += pstate * bio.getRev()[curr][lact][mim][0][sea+w];
									//	mkv.getRepl_open_iof()[curr][mim][t] += pstate * bio.getiofc[curr][lact][mim][0][sea+w];
									//	mkv.getRepl_open_rpo()[curr][mim][t] += pstate * rpo.getCow()[curr][lact][mim][0][sea+w];
									//	mkv.getRepl_open_mlk()[curr][mim][t] += pstate * bio.getMilkyield()[curr][lact][mim][0][sea];
									//	mkv.getRepl_open_cnt()[curr][mim][t] += pstate;
	
									//	replace_file(t,curr,lact,mim,mpreg,pstate,replfrac);
	
	//System.out.println("\n[1] t:%d curr:%d lact:%d mim:%d rpo:%.2f pstate:%f",
	//	  t,curr,lact,mim,rpo.getCow()[curr][lact][mim][0][sea+w],pstate);
	//if(t == 1 && mim <= 11) stop(0);
									}
								}
								else {System.out.println("\n decision error 1"); Utility.stop(1); }
							}
							else if(mim >= bio.getVwp()[lact3] && mim <= bio.getLbm())
							{
	
								p1 = bio.getHeatdetectionrate()[mim][sea] * bio.getConceptionrate()[lact][mim][sea];
								p2 = bio.getHeatdetectionrate()[mim][sea];
								pheat1 = Utility.min2(p2,1.0);								//probility of 1st service
								pheat = Utility.min2(p2 + (1 - p1) * p2 * 9.4 / 21,1.0);		//probility of service in a month
								ppreg = p1 + (1 - p1) * p1 * 9.4 / 21;				//probability of pregnancy in a month
								delig = 1 + (1 - p1) * 9.4 / 21;					//eligible ovulations in a month
								breedcost = pheat * (bio.getP_semen() + bio.getTimebreeding() / 60 * bio.getP_labor());
								reprocost = bio.getReprocost()[mim][sea];				//cow/month
								if(decision == 2)									//decision: serve
								{													//  at end of stage
									mkv.getKeep()[curr][t] += pstate;
									for(next = lo; next <= hi; next++)				//probability in next stage
									{
										mkv.getCow()[next][lact][mim+1][0][jj] += 
											bio.getPnext()[curr][lact][mim][0][next] * pstate * (1 - pid) * (1 - ppreg);
										mkv.getCow()[next][lact][mim+1][1][jj] += 
											bio.getPnext()[curr][lact][mim][0][next] * pstate * (1 - pid) * ppreg;
										mkv.getFs()[next][lact][mim+1][jj] += 
											bio.getPnext()[curr][lact][mim][0][next] * pfs * (1 - pid) * (1 - pheat);
									}
	
	//if(pfs == 0 || bio.getheatdetectionrate[mim][sea] == 0)// || pstate < pfs || mkv.getCow()[next][lact][mim+1][0][jj] < mkv.getFs()[next][lact][mim+1][jj])
	//{
	//System.out.println("\nMARKOV_CHAIN*curr:%d lact:%d mim:%d sea:%d pstate:%.15f pfs:%.15f heat:%f  ppreg:%f",curr,lact,mim,sea,pstate,pfs,bio.getheatdetectionrate[mim][sea]*Constants.STAGE/21,ppreg);
	//stop(0);  
	//}
	
									if(all == 1)
									{
										mkv.getOpen()[curr][t] += pstate;
										mkv.getPind()[curr][t] += pid * pstate;
										mkv.getCash()[t] += pstate * (bio.getRev()[curr][lact][mim][0][sea+w] * df0
													 - reprocost * df0
											         - (1 - pid) * breedcost * df1			//service just after involuntary culling
												     + pid * (- p_loss + p_sell1) * df1);
										mkv.getRev()[t] += pstate * bio.getRev()[curr][lact][mim][0][sea+w] * df0;
										mkv.getMilkyield()[mim][t] += pstate * bio.getMilkyield()[curr][lact][mim][0][sea];
										mkv.getMilkyield_cnt()[mim][t] += pstate;
										if(dry == 1) mkv.getFeedcost()[0][t] += pstate * bio.getDmi()[curr][lact][mim][0][sea] * bio.getP_tmrdmdry() * df0;
										if(dry == 0) mkv.getFeedcost()[1][t] += pstate * bio.getDmi()[curr][lact][mim][0][sea] * bio.getP_tmrdmlac() * df0;
										if(dry == 0) mkv.getDmi()[mim][t] += pstate * bio.getDmi()[curr][lact][mim][0][sea];
										if(dry == 0) mkv.getDmi_cnt()[mim][t] += pstate;
										if(dry == 0) mkv.getMilking()[t] += pstate;
										mkv.getCowsales()[t] += pstate * pid * p_sell1 * df1;	//cows sold at end of stage
										mkv.getLoss()[t] += pstate * pid * p_loss * df1;		//loss due to involuntary culling
										mkv.getLact()[lact][t] += pstate;
										mkv.getDim()[t] += pstate * ((double)mim * Constants.STAGE - (Constants.STAGE / 2));
	
									//	mkv.getLaborcost()[t] += pstate * bio.getlaborcostcows[mim][0] * df0;	//fresh,lact,dry cows
									//	mkv.getLaborcost()[t] += pstate * bio.getlaborcostheat[mim] * df0;	//heat detection or fertility drugs
									//	mkv.getLaborcost()[t] += pstate * (1 - pid) * pheat * bio.getlaborcostbreeding; //labor
	
										mkv.getLaborcost()[t] += pstate * bio.getTimecows()[curr][lact][mim][0][sea] / 60 * bio.getP_labor() * df0;
										mkv.getLaborcost()[t] += pstate * (1 - pid) * pheat * bio.getTimebreeding() / 60 * bio.getP_labor() * df1;
										mkv.getLabortime()[t] += pstate * bio.getTimecows()[curr][lact][mim][0][sea];
										mkv.getLabortime()[t] += pstate * (1 - pid) * pheat * bio.getTimebreeding();
										mkv.getBreedinglaborcost()[t] += pstate * (1 - pid) * pheat * bio.getTimebreeding() / 60 * bio.getP_labor() * df1;	//labor
										mkv.getBreedingsupplycost()[t] += pstate * (reprocost + (1 - pid) * pheat * bio.getP_semen()) * df1;			//drugs + semen
										mkv.getVetcost()[t] += pstate * bio.getVetcost()[mim] * df0;
										mkv.getEligible_serv()[mim][t] += pstate * (1 - pid) * delig;
										mkv.getEligible_bred()[mim][t] += pstate * (1 - pid) * pheat;
										mkv.getEligible_preg()[mim][t] += pstate * (1 - pid) * ppreg;
										mkv.getNewpreg()[curr][lact][t] += pstate * (1 - pid) * ppreg;
	
	//if(lact == 1)
	//System.out.println("\nMARKOV_CHAIN*curr:%d lact:%d mim:%d decision:%d pstate:%.15f ppreg:%f",curr,lact,mim,decision,pstate,ppreg);
	//if(lact == 1) stop(0);
	
	
										mkv.getDaystoconc()[curr][lact][mim][t] += pstate * (1 - pid) * ppreg;
										mkv.getPind_mim()[mim][t] += pstate * pid;
										mkv.getStatus()[mim][0] += pstate;
										mkv.getDysopen()[t] += pstate * (double)mim;
										mkv.getOthercost()[t] += pstate * bio.getP_other_lact() * Constants.STAGE * df0;
										//pregnancy value: at start of state, also for mim == bio.getlbm+1
										if(mim >= bio.getVwp()[lact3] + 1) mkv.getPregval()[curr][lact][mim][0][t] += pstate * rpo.getPrg()[curr][lact][mim][0][sea+w];
										if(mim >= bio.getVwp()[lact3] + 1) mkv.getPregval_cnt()[curr][lact][mim][0][t] += pstate;
										//insemination value: serv - keep in this state
										mkv.getInsval()[mim][t] += pstate * (1- pid) * ins.getCow()[curr][lact][mim][sea+w];
										if(mim == bio.getVwp()[lact3]) mkv.getBreedings1()[t] += pstate * (1 - pid);
										mkv.getHz_open()[mim][0] += pstate;				//enter state
										mkv.getHz_open()[mim][2] += pstate * pid;		//pind open: at end of month
										mkv.getDvwp()[curr][lact][mim][t] += pvw;		//days to start breeding period:									
										mkv.getDfs()[curr][lact][mim][t] += pfs * (1 - pid) * pheat;//days to first service: cows that abort become open but are not eligible
									
	
	//System.out.println("\nt:%d curr:%d lact:%d mim:%d pfs:%f pheat:%f  dfs[]:%f  dfs_cnt[]:%f  dfs:%.1f %.1f",
	//	   	   t,curr,lact,mim,pfs,pheat,mkv.getDfs()[curr][t],mkv.getDfs_cnt()[curr][t],mkv.getDfs()[curr][t]/mkv.getDfs_cnt()[curr][t],(double)mim * Constants.STAGE);
	
	
										mkv.getDlbm()[curr][lact][mim][t] += pstate * (1 - pid);//days to last service
									
										//cost of an extra day open
										mkv.getDysopencost()[curr][mim][t] += pstate * del.getOpen()[curr][lact][mim][mim+1][sea] / Constants.STAGE;
										mkv.getDysopencost_cnt()[curr][mim][t] += pstate;
										//cost of an extra day later pregnant (Monsanto question)
										if(mim > bio.getVwp()[lact3]) mkv.getDyspregcost()[curr][mim-1][t] += pstate * del.getPreg()[curr][lact][mim][mim+1][sea] / Constants.STAGE;
										if(mim > bio.getVwp()[lact3]) mkv.getDyspregcost_cnt()[curr][mim-1][t] += pstate;
	
	//if(mim > 12)
	//{
	//System.out.println("\nt:%d [serv] curr:%d lact:%d mim:%d  open:%.2f  preg:%.2f  dysopencost:%.3f  dyspregcost:%.3f",
	//	  t,curr,lact,mim,del.getOpen()[curr][lact][mim][mim+1][sea]/Constants.STAGE,del.getPreg()[curr][lact][mim][mim+1][sea]/Constants.STAGE,
	//	  mkv.getDysopencost()[curr][mim][t]/mkv.getDysopencost_cnt()[curr][mim][t],
	//	  mkv.getDyspregcost()[curr][mim-1][t]/mkv.getDyspregcost_cnt()[curr][mim-1][t]
	//	  );
	
	//if(del.getPreg()[curr][lact][mim][mim+1][sea]/Constants.STAGE < -5) stop(0);
	//}
										if(rpo.getCow()[curr][lact][mim][0][sea+w] != -Constants.X) mkv.getPind_rpo()[curr][lact][mim][0] += pstate * pid * rpo.getCow()[curr][lact][mim][0][sea+w];
										if(rpo.getCow()[curr][lact][mim][0][sea+w] != -Constants.X) mkv.getPind_rpo_cnt()[curr][lact][mim][0] += pstate * pid;
	
	//if(all == 1)
	//{
	//System.out.println("\n %d %d %d %d :%f  %f  %f",curr,lact,mim,sea+w,rpo.getCow()[curr][lact][mim][0][sea+w],mkv.getPind_rpo()[curr][lact][mim][0],mkv.getPind_rpo_cnt()[curr][lact][mim][0]);
	//if(rpo.getCow()[curr][lact][mim][0][sea+w] == -XXXX) stop(1);
	//}
	
										//do we ever keep open cows that do not pay for themselves?
										if(bio.getRev()[curr][lact][mim][0][sea+w] < 0.0) mkv.getKeep_rev0()[t] += pstate;
									}
								}
								else if(decision == 1)							//decision: keep, do not serve
								{
									mkv.getKeep()[curr][t] += pstate;
									for(next = lo; next <= hi; next++)
									{
										mkv.getCow()[next][lact][mim+1][0][jj] +=
											bio.getPnext()[curr][lact][mim][0][next] * pstate * (1 - pid);
										mkv.getFs()[next][lact][mim+1][jj] += 
											bio.getPnext()[curr][lact][mim][0][next] * pfs * (1 - pid);
										mkv.getVw()[next][lact][mim+1][jj] += 
											bio.getPnext()[curr][lact][mim][0][next] * pvw * (1 - pid);
									}
									if(all == 1)
									{
										mkv.getOpen()[curr][t] += pstate;
										mkv.getPind()[curr][t] += pstate * pid;
										mkv.getCash()[t] += pstate * (bio.getRev()[curr][lact][mim][0][sea+w] * df0
												  	 + pid * (- p_loss + p_sell1) * df1);
										mkv.getRev()[t] += pstate * bio.getRev()[curr][lact][mim][0][sea+w] * df0;
										mkv.getMilkyield()[mim][t] += pstate * bio.getMilkyield()[curr][lact][mim][0][sea];
										mkv.getMilkyield_cnt()[mim][t] += pstate;
										if(dry == 1) mkv.getFeedcost()[0][t] += pstate * bio.getDmi()[curr][lact][mim][0][sea] * bio.getP_tmrdmdry() * df0;
										if(dry == 0) mkv.getFeedcost()[1][t] += pstate * bio.getDmi()[curr][lact][mim][0][sea] * bio.getP_tmrdmlac() * df0;
										if(dry == 0) mkv.getDmi()[mim][t] += pstate * bio.getDmi()[curr][lact][mim][0][sea];
										if(dry == 0) mkv.getDmi_cnt()[mim][t] += pstate;
										if(dry == 0) mkv.getMilking()[t] += pstate;
										mkv.getCowsales()[t] += pstate * pid * p_sell1 * df1;	//cows sold at end of stage
										mkv.getLoss()[t] += pstate * pid * p_loss * df1;		//loss due to involuntary culling
										mkv.getLact()[lact][t] += pstate;
										mkv.getDim()[t] += pstate * ((double)mim * Constants.STAGE - (Constants.STAGE / 2));
										mkv.getEligible_keep()[mim][t] += pstate * (1 - pid) * delig;
										mkv.getLaborcost()[t] += pstate * bio.getTimecows()[curr][lact][mim][0][sea] / 60 * bio.getP_labor() * df0;
										mkv.getLabortime()[t] += pstate * bio.getTimecows()[curr][lact][mim][0][sea];
										mkv.getVetcost()[t] += pstate * bio.getVetcost()[mim] * df0;
										mkv.getPind_mim()[mim][t] += pstate * pid;
										mkv.getStatus()[mim][0] += pstate;
										mkv.getDysopen()[t] += pstate * (double)(mim);
										mkv.getOthercost()[t] += pstate * bio.getP_other_lact() * Constants.STAGE * df0;
										//insemination value: serv - keep in this state
										mkv.getInsval()[mim][t] += pstate * (1 - pid) * ins.getCow()[curr][lact][mim][sea+w];
										mkv.getHz_open()[mim][0] += pstate;				//enter state
										mkv.getHz_open()[mim][2] += pstate * pid;		//pind open: at end of month
										//cost of an extra day later pregnant
										if(mim > bio.getVwp()[lact3]) mkv.getDyspregcost()[curr][mim-1][t] += pstate * del.getPreg()[curr][lact][mim][mim+1][sea] / Constants.STAGE;
										if(mim > bio.getVwp()[lact3]) mkv.getDyspregcost_cnt()[curr][mim-1][t] += pstate;
										//cost of an extra day open
										mkv.getDysopencost()[curr][mim][t] += pstate * del.getOpen()[curr][lact][mim][mim+1][sea] / Constants.STAGE;
										mkv.getDysopencost_cnt()[curr][mim][t] += pstate;
										if(rpo.getCow()[curr][lact][mim][0][sea+w] != -Constants.X) mkv.getPind_rpo()[curr][lact][mim][0] += pstate * pid * rpo.getCow()[curr][lact][mim][0][sea+w];
										if(rpo.getCow()[curr][lact][mim][0][sea+w] != -Constants.X) mkv.getPind_rpo_cnt()[curr][lact][mim][0] += pstate * pid;
	
	//if(mim > bio.getvwp) 
	//{
	//System.out.println("\nt:%d [keep] curr:%d lact:%d mim:%d  prg:%.2f  dysopencost:%.3f",
	//	  t,curr,lact,mim,del.getPreg()[curr][lact][mim][mim+1][sea]/Constants.STAGE,mkv.getDysopencost()[curr][mim-1][t]/mkv.getDysopencost_cnt()[curr][mim-1][t]);
	//if(del.getPreg()[curr][lact][mim][mim+1][sea]/Constants.STAGE < -5) stop(0);
	//}
										//do we ever keep open cows that do not pay for themselves?
										if(bio.getRev()[curr][lact][mim][0][sea+w] < 0.0) mkv.getKeep_rev0()[t] += pstate;
									}
								}
								else if(decision == 0)							//decision: replace
								{
									if(all == 1)
									{
										mkv.getRepl()[curr][t] += pstate;
										mkv.getRepl_mim()[mim][t] += pstate;
										mkv.getCash()[t] += pstate * p_sell0 * df0;
										mkv.getCowsales()[t] += pstate * p_sell0 * df0;	//cows sold at begin of stage
										mkv.getHz_open()[mim][0] += pstate;				//enter state
										mkv.getHz_open()[mim][1] += pstate;				//repl open: at start of month
	
										mkv.getRepl_open_rev()[curr][mim][t] += pstate * bio.getRev()[curr][lact][mim][0][sea+w];
										mkv.getRepl_open_iof()[curr][mim][t] += pstate * bio.getIofc()[curr][lact][mim][0][sea+w];
										mkv.getRepl_open_rpo()[curr][mim][t] += pstate * rpo.getCow()[curr][lact][mim][0][sea+w];
										mkv.getRepl_open_mlk()[curr][mim][t] += pstate * bio.getMilkyield()[curr][lact][mim][0][sea];
										mkv.getRepl_open_cnt()[curr][mim][t] += pstate;
									//	replace_file(t,curr,lact,mim,mpreg,pstate,replfrac);
	
	
	//if(rpo.getCow()[curr][lact][mim][0][sea+w] > 0)
	//{
	//System.out.println("\n[2] t:%d curr:%d lact:%d mim:%d rpo:%.2f pstate:%f pol:%d",
	//	  t,curr,lact,mim,rpo.getCow()[curr][lact][mim][0][sea+w],pstate,pol.getCow()[curr][lact][mim][0][sea+w]);
	//if(t == 1 && mim <= 11) stop(0);
	//stop(0);
	//}
	
									}
								}
								else {System.out.println("\n decision error 2"); Utility.stop(1); }
							}
							else if(mim <= bio.getVwp()[lact3] - 1 || mim >= bio.getLbm() + 1)
							{
								if(lact == 1 && mim == 1)
								{
									//entering heifers calculated at end of stage
								}
								else if(decision == 1)							//decision: keep
								{
									mkv.getKeep()[curr][t] += pstate;
									for(next = lo; next <= hi; next++)
									{
										mkv.getCow()[next][lact][mim+1][0][jj] += 
											bio.getPnext()[curr][lact][mim][0][next] * pstate * (1 - pid);
										mkv.getFs()[next][lact][mim+1][jj] += 
											bio.getPnext()[curr][lact][mim][0][next] * pfs * (1 - pid);
									}
									if(all == 1)
									{
										mkv.getOpen()[curr][t] += pstate;
										mkv.getPind()[curr][t] += pid * pstate;
										mkv.getCash()[t] += pstate * (bio.getRev()[curr][lact][mim][0][sea+w] * df0
												  	 + pid * (- p_loss + p_sell1) * df1);
										mkv.getRev()[t] += pstate * bio.getRev()[curr][lact][mim][0][sea+w] * df0;
										mkv.getMilkyield()[mim][t] += pstate * bio.getMilkyield()[curr][lact][mim][0][sea];
										mkv.getMilkyield_cnt()[mim][t] += pstate;
										if(dry == 1) mkv.getFeedcost()[0][t] += pstate * bio.getDmi()[curr][lact][mim][0][sea] * bio.getP_tmrdmdry() * df0;
										if(dry == 0) mkv.getFeedcost()[1][t] += pstate * bio.getDmi()[curr][lact][mim][0][sea] * bio.getP_tmrdmlac() * df0;
										if(dry == 0) mkv.getDmi()[mim][t] += pstate * bio.getDmi()[curr][lact][mim][0][sea];
										if(dry == 0) mkv.getDmi_cnt()[mim][t] += pstate;
										if(dry == 0) mkv.getMilking()[t] += pstate;
										mkv.getCowsales()[t] += pstate * pid * p_sell1 * df1;	//cows sold at end of stage
										mkv.getLoss()[t] += pstate * pid * p_loss * df1;		//loss due to involuntary culling
										mkv.getLact()[lact][t] += pstate;
										mkv.getDim()[t] += pstate * ((double)mim * Constants.STAGE - (Constants.STAGE / 2));
										mkv.getLaborcost()[t] += pstate * bio.getTimecows()[curr][lact][mim][0][sea] / 60 * bio.getP_labor() * df0;
										mkv.getLabortime()[t] += pstate * bio.getTimecows()[curr][lact][mim][0][sea];
										mkv.getVetcost()[t] += pstate * bio.getVetcost()[mim] * df0;
										if(mim == 1) mkv.getCalves_born()[t] += pstate;
										if(mim == 1) mkv.getCash()[t] += pstate * bio.getP_calf() * df0;
										if(mim == 1) mkv.getCalved()[t] += pstate;
										mkv.getPind_mim()[mim][t] += pstate * pid;
										mkv.getStatus()[mim][0] += pstate;
										mkv.getDysopen()[t] += pstate * (double)(mim);
										mkv.getOthercost()[t] += pstate * bio.getP_other_lact() * Constants.STAGE * df0;
										mkv.getHz_open()[mim][0] += pstate;				//enter state
										mkv.getHz_open()[mim][2] += pstate * pid;		//pind open: at end of month
										//pregnancy value open cow: at start of state, see DP_POLICY
										if(mim == bio.getLbm() + 1) mkv.getPregval()[curr][lact][mim][0][t] += pstate * rpo.getPrg()[curr][lact][mim][0][sea+w];
										if(mim == bio.getLbm() + 1) mkv.getPregval_cnt()[curr][lact][mim][0][t] += pstate;
										if(rpo.getCow()[curr][lact][mim][0][sea+w] != -Constants.X) mkv.getPind_rpo()[curr][lact][mim][0] += pstate * pid * rpo.getCow()[curr][lact][mim][0][sea+w];
										if(rpo.getCow()[curr][lact][mim][0][sea+w] != -Constants.X) mkv.getPind_rpo_cnt()[curr][lact][mim][0] += pstate * pid;
										//do we ever keep open cows that do not pay for themselves?
										if(bio.getRev()[curr][lact][mim][0][sea+w] < 0.0) mkv.getKeep_rev0()[t] += pstate;
									}
								}
								else if(decision == 0)							//decision: replace
								{
									if(all == 1)
									{
	
										mkv.getRepl()[curr][t] += pstate;
										mkv.getCash()[t] += pstate * p_sell0 * df0;
										mkv.getCowsales()[t] += pstate * p_sell0 * df0;	//cows sold at begin of stage
										if(mim == 1) mkv.getCalves_born()[t] += pstate;
										if(mim == 1) mkv.getCash()[t] += pstate * bio.getP_calf() * df0;
										if(mim == 1) mkv.getCalved()[t] += pstate;
										mkv.getRepl_mim()[mim][t] += pstate;
										mkv.getHz_open()[mim][0] += pstate;				//enter state
										mkv.getHz_open()[mim][1] += pstate;				//repl open: at start of month
	
										mkv.getRepl_open_rev()[curr][mim][t] += pstate * bio.getRev()[curr][lact][mim][0][sea+w];
										mkv.getRepl_open_iof()[curr][mim][t] += pstate * bio.getIofc()[curr][lact][mim][0][sea+w];
										mkv.getRepl_open_rpo()[curr][mim][t] += pstate * rpo.getCow()[curr][lact][mim][0][sea+w];
										mkv.getRepl_open_mlk()[curr][mim][t] += pstate * bio.getMilkyield()[curr][lact][mim][0][sea];
										mkv.getRepl_open_cnt()[curr][mim][t] += pstate;
	
									//	replace_file(t,curr,lact,mim,mpreg,pstate,replfrac);
	
	
	//if(rpo.getCow()[curr][lact][mim][0][sea+w] > 0)
	//{
	//System.out.println("\n[1] t:%d curr:%d lact:%d mim:%d rpo:%.2f pstate:%f pol:%d",
	//	  t,curr,lact,mim,rpo.getCow()[curr][lact][mim][0][sea+w],pstate,pol.getCow()[curr][lact][mim][0][sea+w]);
	//if(t == 1 && mim <= 11) stop(0);
	//stop(0);
	//}
									}
								}
								else 
								{
									System.out.println("\n decision error 3: decision:"+decision); 
									System.out.println("\n[2] t:"+t+" curr:"+curr+" lact:"+lact+" mim:"+mim+" rpo:"+rpo.getCow()[curr][lact][mim][0][sea+w]+" pstate:"+pstate);
									Utility.stop(1); 
								}
							}
						}//end if(pstate > 0.0)
	//*** PREGNANT COWS **********************************************************************
						for(mpreg = 1; mpreg <= Constants.MPREG; mpreg++) if(mim - mpreg >= bio.getVwp()[lact3] && mim - mpreg <= bio.getLbm())
						{
							pstate = mkv.getCow()[curr][lact][mim][mpreg][ii];
							if(pstate > 0.00000000000001 && alldelayed == 0)
							{
								decision = pol.getCow()[curr][lact][mim][mpreg][sea+w];
								if(decision == Constants.POL) { System.out.println("\n decision error preg: "+curr+" "+lact+" "+mim+" "+(mpreg+sea+w)+" : "+pol.getCow()[curr][lact][mim][mpreg][sea+w]); Utility.stop(1); }
								if(decision != 0 && decision != 1) { System.out.println("\n decision error preg ="+decision); Utility.stop(1); }
	
	
	//System.out.println("\n[Constants.MPREG] t:%d curr:%d lact:%d mim:%d rpo:%.2f pstate:%f pol:%d",
	//	  t,curr,lact,mim,rpo.getCow()[curr][lact][mim][mpreg][sea+w],pstate,pol.getCow()[curr][lact][mim][mpreg][sea+w]);
	//if(rpo.getCow()[curr][lact][mim][mpreg][sea+w] > 0 && decision == 0)
	//{
	//System.out.println("\n[0] t:%d curr:%d lact:%d mim:%d rpo:%.2f pstate:%f pol:%d",
	//	  t,curr,lact,mim,rpo.getCow()[curr][lact][mim][0][sea+w],pstate,pol.getCow()[curr][lact][mim][0][sea+w]);
	//stop(0);
	//}
	
	
								p_sell0 = bio.getSel()[lact][mim][sea+w];
								if(sea == Constants.SEA) p_sell2 = bio.getSel()[lact][mim][1+w];//salvage value end of stage
								if(sea <  Constants.SEA) p_sell2 = bio.getSel()[lact][mim][sea+1+w];
								p_sell1 = p_sell2 * bio.getInvcullpct();	
	
								if(bio.getMilkyield()[curr][lact][mim][mpreg][sea] > 0.01) dry = 0; else dry = 1;
	
								pid = bio.getPind()[lact][mim][sea];
								abort = bio.getAbort()[mpreg];
								if(all == 1)
								{
									if(rpo.getCow()[curr][lact][mim][mpreg][sea+w] != -Constants.X) mkv.getRpo()[t] += pstate * rpo.getCow()[curr][lact][mim][mpreg][sea+w];
									if(rpo.getCow()[curr][lact][mim][mpreg][sea+w] != -Constants.X) mkv.getRpo_cnt()[t] += pstate;
									if(rpo.getCow()[curr][lact][mim][mpreg][sea+w] > 0) mkv.getValue()[t] += pstate * (rpo.getCow()[curr][lact][mim][mpreg][sea+w] + p_sell0);
									else mkv.getValue()[t] += pstate * p_sell0;
									mkv.getValue_cnt()[t] += pstate;
								}
								if(mpreg <= Constants.MPREG - 1)
								{
									if(decision == 1)							//decision: keep
									{
										mkv.getKeep()[curr][t] += pstate;
										lo = bio.getLo()[curr][lact][mim][0];
										hi = bio.getHi()[curr][lact][mim][0];
										for(next = lo; next <= hi; next++)
										{
											mkv.getCow()[next][lact][mim+1][mpreg+1][jj] += 
												bio.getPnext()[curr][lact][mim][0][next] * pstate * (1 - pid) * (1 - abort);
											mkv.getCow()[next][lact][mim+1][0][jj] += 
												bio.getPnext()[curr][lact][mim][0][next] * pstate * (1 - pid) * abort;
										}
	
										if(all == 1)
										{
											mkv.getPreg()[curr][t] += pstate;
											mkv.getPind()[curr][t] += pid * pstate;
											//Monte Carlo try-out: 	if(mc == 1) mkv.getCow()[curr][lact][mim+1][mpreg+1][jj] += ran() <= (1 - pid) ? 1 : 0;
											mkv.getCash()[t] += pstate * (bio.getRev()[curr][lact][mim][mpreg][sea+w] * df0
												 		+ pid * (- p_loss + p_sell1) * df1);
											mkv.getRev()[t] += pstate * bio.getRev()[curr][lact][mim][mpreg][sea+w] * df0;
											mkv.getMilkyield()[mim][t] += pstate * bio.getMilkyield()[curr][lact][mim][mpreg][sea];
											mkv.getMilkyield_cnt()[mim][t] += pstate;
											if(dry == 1) mkv.getFeedcost()[0][t] += pstate * bio.getDmi()[curr][lact][mim][mpreg][sea] * bio.getP_tmrdmdry() * df0;
											if(dry == 0) mkv.getFeedcost()[1][t] += pstate * bio.getDmi()[curr][lact][mim][mpreg][sea] * bio.getP_tmrdmlac() * df0;
											if(dry == 0) mkv.getDmi()[mim][t] += pstate * bio.getDmi()[curr][lact][mim][mpreg][sea];
											if(dry == 0) mkv.getDmi_cnt()[mim][t] += pstate;
											if(dry == 0) mkv.getMilking()[t] += pstate;
											mkv.getAbortions()[t] += pstate * abort;
											mkv.getAbortions_cnt()[t] += pstate;
											mkv.getLact()[lact][t] += pstate;
											mkv.getDim()[t] += pstate * ((double)mim * Constants.STAGE - (Constants.STAGE / 2));
											mkv.getLaborcost()[t] += pstate * bio.getTimecows()[curr][lact][mim][mpreg][sea] / 60 * bio.getP_labor() * df0;
											mkv.getLabortime()[t] += pstate * bio.getTimecows()[curr][lact][mim][mpreg][sea];
											mkv.getVetcost()[t] += pstate * bio.getVetcost()[mim] * df0;
											mkv.getCowsales()[t] += pstate * pid * p_sell1 * df1;	//cows sold at end of stage
											mkv.getLoss()[t] += pstate * pid * p_loss * df1;		//loss due to involuntary culling
											mkv.getPind_mim()[mim][t] += pstate * pid;
											mkv.getStatus()[mim][1] += pstate;
											if(dry == 0) mkv.getOthercost()[t] += pstate * bio.getP_other_lact() * Constants.STAGE * df0;
											else mkv.getOthercost()[t] += pstate * bio.getP_other_dry() * Constants.STAGE * df0;
											//pregnancy value: extra pay-off because this cow is pregnant vs open in this stage
											mkv.getPregval()[curr][lact][mim][mpreg][t] += pstate * rpo.getPrg()[curr][lact][mim][mpreg][sea+w];
											mkv.getPregval_cnt()[curr][lact][mim][mpreg][t] += pstate;
											mkv.getHz_preg()[mim][0] += pstate;				//enter state
											mkv.getHz_preg()[mim][2] += pstate * pid;		//pind preg: at end of month
											if(rpo.getCow()[curr][lact][mim][mpreg][sea+w] != -Constants.X) mkv.getPind_rpo()[curr][lact][mim][mpreg] += pstate * pid * rpo.getCow()[curr][lact][mim][mpreg][sea+w];
											if(rpo.getCow()[curr][lact][mim][mpreg][sea+w] != -Constants.X) mkv.getPind_rpo_cnt()[curr][lact][mim][mpreg] += pstate * pid;
	
											if(bio.getMilkyield()[curr][lact][mim][mpreg][sea] > 0.01 && bio.getMilkyield()[curr][lact][mim+1][mpreg+1][sea1] <= 0.01)
											{	//dry-off yield
												mkv.getDryoffyield()[t] += pstate * bio.getMilkyield()[curr][lact][mim][mpreg][sea] / Constants.STAGE;
												mkv.getDryoffyield_cnt()[t] += pstate;
											//	if(t == 1)
											//		System.out.println("\n t:%d curr:%d lact:%d mim:%d mpreg:%d yield:%f",t,curr,lact,mim,mpreg,bio.getMilkyield()[curr][lact][mim][mpreg][sea]/Constants.STAGE);
											}
										}
									}
									else if(decision == 0)							//decision: replace
									{
										if(all == 1)
										{
											mkv.getRepl()[curr][t] += pstate;
											mkv.getCash()[t] += pstate * p_sell0 * df0;
											mkv.getCowsales()[t] += pstate * p_sell0 * df0;	//cows sold at begin of stage
											mkv.getRepl_mim()[mim][t] += pstate;
											mkv.getHz_preg()[mim][0] += pstate;				//enter state
											mkv.getHz_preg()[mim][1] += pstate;				//repl preg: at start of month
	
	//System.out.println("\n[1] t:%d curr:%d lact:%d mim:%d rpo:%.2f pstate:%.10f pol:%d decision:%d",
	//	  t,curr,lact,mim,rpo.getCow()[curr][lact][mim][0][sea+w],pstate,pol.getCow()[curr][lact][mim][mpreg][sea+w],decision);
	//stop(0);
	
										//	replace_file(t,curr,lact,mim,mpreg,pstate,replfrac);
										}
									}
									else {System.out.println("\n decision error 4"); Utility.stop(1); }
								}
								else if(mpreg == Constants.MPREG && lact <= bio.getMaxlact() - 1)		//end of lactation
								{
									if(decision == 1)							//decision: keep
									{
										mkv.getKeep()[curr][t] += pstate;
										lo = bio.getLo()[curr][lact][mim][sea];
										hi = bio.getHi()[curr][lact][mim][sea];
										for(next = lo; next <= hi; next++)
										{
											mkv.getCow()[next][lact+1][1][0][jj] += 
												bio.getPnext()[curr][lact][mim][sea][next] * pstate * (1 - pid);
											mkv.getFs()[next][lact+1][1][jj] += 
												bio.getPnext()[curr][lact][mim][sea][next] * pstate * (1 - pid);
										}
										if(all == 1)
										{
											mkv.getPreg()[curr][t] += pstate;
											mkv.getPind()[curr][t] += pid * pstate;
											mkv.getCash()[t] += pstate * (bio.getRev()[curr][lact][mim][mpreg][sea+w] * df0
													 	 + pid * (- p_loss + p_sell1) * df1);
											mkv.getRev()[t] += pstate * bio.getRev()[curr][lact][mim][mpreg][sea+w] * df0;
											mkv.getMilkyield()[mim][t] += pstate * bio.getMilkyield()[curr][lact][mim][mpreg][sea];
											mkv.getMilkyield_cnt()[mim][t] += pstate;
											if(dry == 1) mkv.getFeedcost()[0][t] += pstate * bio.getDmi()[curr][lact][mim][mpreg][sea] * bio.getP_tmrdmdry() * df0;
											if(dry == 0) mkv.getFeedcost()[1][t] += pstate * bio.getDmi()[curr][lact][mim][mpreg][sea] * bio.getP_tmrdmlac() * df0;
											if(dry == 0) mkv.getDmi()[mim][t] += pstate * bio.getDmi()[curr][lact][mim][mpreg][sea];
											if(dry == 0) mkv.getDmi_cnt()[mim][t] += pstate;
											if(dry == 0) mkv.getMilking()[t] += pstate;
											mkv.getCowsales()[t] += pstate * pid * p_sell1 * df1;	//all cows sold at end of stage
											mkv.getLoss()[t] += pstate * pid * p_loss * df1;		//loss due to involuntary culling
											mkv.getClvinterval()[lact][mim][t] += pstate * (1 - pid);		//calving should be at start of next t, but otherwise t=1 does not have a calving,
											mkv.getLact()[lact][t] += pstate;						//		because mim in previous lactation is not remembered
											mkv.getDim()[t] += pstate * ((double)mim * Constants.STAGE - (Constants.STAGE / 2));
											mkv.getLaborcost()[t] += pstate * bio.getTimecows()[curr][lact][mim][mpreg][sea] / 60 * bio.getP_labor() * df0;
											mkv.getLabortime()[t] += pstate * bio.getTimecows()[curr][lact][mim][mpreg][sea];
											mkv.getVetcost()[t] += pstate * bio.getVetcost()[mim] * df0;
											mkv.getPind_mim()[mim][t] += pstate * pid;
											mkv.getStatus()[mim][1] += pstate;
											mkv.getOthercost()[t] += pstate * bio.getP_other_dry() * Constants.STAGE * df0;
											//pregnancy value: extra pay-off because this cow is pregnant vs open in this stage
											mkv.getPregval()[curr][lact][mim][mpreg][t] += pstate * rpo.getPrg()[curr][lact][mim][mpreg][sea+w];
											mkv.getPregval_cnt()[curr][lact][mim][mpreg][t] += pstate;
											mkv.getHz_preg()[mim][0] += pstate;				//enter state
											mkv.getHz_preg()[mim][2] += pstate * pid;		//pind preg: at end of month
											if(rpo.getCow()[curr][lact][mim][mpreg][sea+w] != -Constants.X) mkv.getPind_rpo()[curr][lact][mim][mpreg] += pstate * pid * rpo.getCow()[curr][lact][mim][mpreg][sea+w];
											if(rpo.getCow()[curr][lact][mim][mpreg][sea+w] != -Constants.X) mkv.getPind_rpo_cnt()[curr][lact][mim][mpreg] += pstate * pid;
										}
									}
									else if(decision == 0)						//decision: replace
									{
										if(all == 1)
										{
											mkv.getRepl()[curr][t] += pstate;
											mkv.getCash()[t] += pstate * p_sell0 * df0;
											mkv.getCowsales()[t] += pstate * p_sell0 * df0;	//cows sold at begin of stage
											mkv.getRepl_mim()[mim][t] += pstate;
											mkv.getHz_preg()[mim][0] += pstate;				//enter state
											mkv.getHz_preg()[mim][1] += pstate;				//repl preg: at start of month
	
	//System.out.println("\n[2] t:%d curr:%d lact:%d mim:%d rpo:%.2f pstate:%.10f pol:%d decision:%d",
	//	  t,curr,lact,mim,rpo.getCow()[curr][lact][mim][0][sea+w],pstate,pol.getCow()[curr][lact][mim][mpreg][sea+w],decision);
	//stop(0);
	
										//	replace_file(t,curr,lact,mim,mpreg,pstate,replfrac);
										}
									}
									else {System.out.println("\n decision error 5"); Utility.stop(1); }
								}
								else if(mpreg == Constants.MPREG && lact == bio.getMaxlact())			//end of lactation Constants.LACT
								{
									decision = 0;	//! Cows going to the bio.getMaxlact()+1 lactation would not be counted as present 
													//  and reduce the fraction present in the herd.
									if(decision == 1)									//decision: keep
									{
										mkv.getKeep()[curr][t] += pstate;
										if(all == 1)
										{
											mkv.getPreg()[curr][t] += pstate;
											mkv.getPind()[curr][t] += pstate * pid;	//all cows culled at end of stage
											mkv.getRepl()[curr][t] += pstate * (1 - pid);	//all cows culled at end of stage
											mkv.getCash()[t] += pstate * (bio.getRev()[curr][lact][mim][mpreg][sea+w] * df0
													  	 + pid * (- p_loss + p_sell1) * df1
														 + (1 - pid) * p_sell2 * df1);
											mkv.getRev()[t] += pstate * bio.getRev()[curr][lact][mim][mpreg][sea+w] * df0;
											mkv.getMilkyield()[mim][t] += pstate * bio.getMilkyield()[curr][lact][mim][mpreg][sea];
											mkv.getMilkyield_cnt()[mim][t] += pstate;
											if(dry == 1) mkv.getFeedcost()[0][t] += pstate * bio.getDmi()[curr][lact][mim][mpreg][sea] * bio.getP_tmrdmdry() * df0;
											if(dry == 0) mkv.getFeedcost()[1][t] += pstate * bio.getDmi()[curr][lact][mim][mpreg][sea] * bio.getP_tmrdmlac() * df0;
											if(dry == 0) mkv.getDmi()[mim][t] += pstate * bio.getDmi()[curr][lact][mim][mpreg][sea];
											if(dry == 0) mkv.getDmi_cnt()[mim][t] += pstate;
											if(dry == 0) mkv.getMilking()[t] += pstate;
											mkv.getCowsales()[t] += pstate * pid * p_sell1 * df1;//all cows sold at end of stage
											mkv.getCowsales()[t] += pstate * (1 - pid) * p_sell2 * df1;
											mkv.getLoss()[t] += pstate * pid * p_loss * df1;	//loss due to involuntary culling
											mkv.getClvinterval()[lact][mim][t] += pstate * (1 - pid);
											mkv.getLact()[lact][t] += pstate;
											mkv.getDim()[t] += pstate * ((double)mim * Constants.STAGE - (Constants.STAGE / 2));
											mkv.getLaborcost()[t] += pstate * bio.getTimecows()[curr][lact][mim][mpreg][sea] / 60 * bio.getP_labor() * df0;
											mkv.getLabortime()[t] += pstate * bio.getTimecows()[curr][lact][mim][mpreg][sea];
											mkv.getVetcost()[t] += pstate * bio.getVetcost()[mim] * df0;
											mkv.getPind_mim()[mim][t] += pstate * pid;
											mkv.getStatus()[mim][1] += pstate;
											mkv.getRepl_mim()[cs0][t] += pstate * (1 - pid);
											mkv.getOthercost()[t] += pstate * bio.getP_other_dry() * Constants.STAGE * df0;
											//pregnancy value: extra pay-off because this cow is pregnant vs open in this stage
											mkv.getPregval()[curr][lact][mim][mpreg][t] += pstate * rpo.getPrg()[curr][lact][mim][mpreg][sea+w];
											mkv.getPregval_cnt()[curr][lact][mim][mpreg][t] += pstate;
											mkv.getHz_preg()[mim][0] += pstate;				//enter state
										//	mkv.getHz_preg()[mim][1] += pstate * (1 - pid);	//repl preg: at end of month
											mkv.getHz_preg()[mim][2] += pstate * pid;		//pind preg: at end of month
											if(rpo.getCow()[curr][lact][mim][mpreg][sea+w] != -Constants.X) mkv.getPind_rpo()[curr][lact][mim][mpreg] += pstate * pid * rpo.getCow()[curr][lact][mim][mpreg][sea+w];
											if(rpo.getCow()[curr][lact][mim][mpreg][sea+w] != -Constants.X) mkv.getPind_rpo_cnt()[curr][lact][mim][mpreg] += pstate * pid;
										}	//see mim > bio.getmaxmim for hz_open
									}
									else if(decision == 0)						//decision: replace
									{
										if(all == 1)
										{
											mkv.getRepl()[curr][t] += pstate;
											mkv.getCash()[t] += pstate * p_sell0 * df0;
											mkv.getCowsales()[t] += pstate * p_sell0 * df0;	//cows sold at begin of stage
											mkv.getRepl_mim()[mim][t] += pstate;
											mkv.getHz_preg()[mim][0] += pstate;				//enter state
											mkv.getHz_preg()[mim][1] += pstate;				//repl preg: at start of month
	
	//System.out.println("\n[3] t:%d curr:%d lact:%d mim:%d rpo:%.2f pstate:%.10f pol:%d decision:%d",
	//	  t,curr,lact,mim,rpo.getCow()[curr][lact][mim][0][sea+w],pstate,pol.getCow()[curr][lact][mim][mpreg][sea+w],decision);
	//stop(0);
	
										//	replace_file(t,curr,lact,mim,mpreg,pstate,replfrac);
										}
									}
									else {System.out.println("\n decision error 6"); Utility.stop(1); }
								}
								else { System.out.println("\n MARKOV_CHAIN pregnancy error"); Utility.stop(1); }
							}//end if(pstate > 0.0)
			//				if(repeat == 2) goto repeat2;	//repeat program from repeat2
						}//end for(mpreg)
					}//end for(mim)
				}//end for(curr)
			}//end for(lact)
	/****** ENTERING HEIFERS *****************************************************/
			//mkv.getheiferpool[t] is the actual fraction heifers brought into the herd:
			//  It is a function of the heifers raised and purchased.
			//  Raised heifers go into heiferpool before heifers are purchased.
			if(all == 1) for(age = Constants.MPREG; age <= Constants.AGE; age++)				//raised heifer calves
				mkv.getHeifers_raised()[t] += mkv.getHeifer()[age][Constants.MPREG][ii];
	
			//LP simulation of starting herd/cow over time
			//	lpt = 0: herd, lpt >= 1: single heifer entering at lpt = t
			if(pool == 0)
			{
			//	if(lpt == t) mkv.getheiferpool[t] = 1.0;				//1.0 heifer entering at t == lpt
			//	else mkv.getheiferpool[t] = 0.0;						//0.0 heifer entering at t != 1
	
				if(t == 1)
				{
					hlp = 0.0; for(curr = 1; curr <= Constants.CURR; curr++) hlp += mkv.getKeep()[curr][t];
					for(curr = 1; curr <= Constants.CURR; curr++) hlp += mkv.getRepl()[curr][t];
					if(hlp == 0.0) mkv.getHeiferpool()[t] = 1.0;		//empty herd .get enter 1 heifer at t = 1
					else mkv.getHeiferpool()[t] = 0.0;					//no heifers entering
				}
				else  mkv.getHeiferpool()[t] = 0.0;
			}
			//simulation: read markov.prn + read heiferpool.txt
	//		else if(pool == 2 || pool == 9)							//2=GA, 9=LP
	//		{
	//			//if(pool == 9 && t > MAXU)
	//			//	if(dpp.getenter[sea] == 1)						//optimal DP decision
	//			//	{
	//			//		hlp = 0.0; for(curr = 1; curr <= Constants.CURR; curr++) hlp += mkv.getKeep()[curr][t];
	//			//		mkv.getheiferpool[t] = 1.0 - hlp;				//fraction heifers entering
	//			//	}
	//			//	else
	//			//		mkv.getheiferpool[t] = 0.0;
	//			System.out.println("\n MARKOV_CHAIN: error entering heifers: GA or LP");
	//			stop(4);
	//		}
			//DP: read markov.prn, calculate heiferpool
			else if(pool == 1)
			{
				if(all == 1) mkv.getHeifers_raised()[t] = 0.0;			//to make heifers_purchased == heiferpool
	
			//	if(bio.getenter[sea+w] == 1)							//########DBC2006: use for LP SME herds
				if(pol.getEnter()[sea+w] == 1)							//normal use
				{
			//		if(t == 1) System.out.println("*bio.getenter*");				//########DBC2006: use for LP SME herds		
					hlp = 0.0; for(curr = 1; curr <= Constants.CURR; curr++) hlp += mkv.getKeep()[curr][t];
					mkv.getHeiferpool()[t] = 1.0 - hlp;					//fraction heifers entering
				}
				else 
			//	if(bio.getenter[sea+w] == 0)							//########DBC2006: use for LP SME herds	
				if(pol.getEnter()[sea+w] == 0)							//normal use
				{
					mkv.getHeiferpool()[t] = 0.0;
				}
				else
				{
					System.out.println("\n MARKOV_CHAIN: pol.getenter: "+pol.getEnter()[sea+w]+" error in ENTERING HEIFERS");
					Utility.stop(1);
				}
			}
			else {System.out.println("\n decision error 7"); Utility.stop(1); }
			mkv.getEnter()[t] = pol.getEnter()[sea+w];
			if(all == 1)
			{
				mkv.getHeifers_purchased()[t] = Utility.max2(mkv.getHeiferpool()[t] - mkv.getHeifers_raised()[t],0.0);
				mkv.getHeifers_sold()[t] = Utility.max2(mkv.getHeifers_raised()[t] - mkv.getHeiferpool()[t],0.0);
				mkv.getCalves_born()[t] += mkv.getHeiferpool()[t];				//calves born to heifers
				mkv.getHeifer()[2][0][jj] = mkv.getCalves_born()[t] * 0.5 * 0.85;	//85% survival rate of 50% calves
				mkv.getFcalfcost()[t] += mkv.getCalves_born()[t] * 0.50 * p_heifer / 24 * df0;//1 month of heifer raising cost
				mkv.getCash()[t] += mkv.getHeiferpool()[t] * - p_heifer * df0;	//DP: purchase all heifers
				mkv.getCash()[t] += mkv.getHeiferpool()[t] * bio.getP_calf() * df0;
			}
			lact = 1; mim = 1; mpreg = 0;
			pid = bio.getPind()[1][1][sea];									//pid at end of stage
			if(sea == Constants.SEA) p_sell1 = bio.getSel()[1][1][1+w] * bio.getInvcullpct(); 
			else p_sell1 = bio.getSel()[1][1][sea+1+w] * bio.getInvcullpct();
			for(curr = 1; curr <= Constants.CURR; curr++)							//Constants.CURR milk yield levels of heifers that enter herd
			{
				pstate = mkv.getHeiferpool()[t] * bio.getPnext()[0][0][0][0][curr];//disperse mkv.getheiferpool[] over Constants.CURR classes
			//	mkv.getCow()[curr][1][1][0][ii] = pstate;					//start of this stage
				lo = bio.getLo()[curr][1][1][0];
				hi = bio.getHi()[curr][1][1][0];
				mkv.getKeep()[curr][t] += pstate;
				for(next = lo; next <= hi; next++)
				{
					mkv.getCow()[next][1][2][0][jj] += 
						bio.getPnext()[curr][1][1][0][next] * pstate * (1 - pid);	//start of next stage, mkvsum next stage < 1.0
				//	mkv.getFs()[next][1][2][jj] = mkv.getCow()[next][1][2][0][jj];
				//	mkv.getVw()[next][1][2][jj] = mkv.getCow()[next][1][2][0][jj];
				}
				if(all == 1)
				{
					mkv.getCash()[t] += pstate * (bio.getRev()[curr][1][1][0][sea+w] * df0	//cost during this stage
							     + pid * (- p_loss + p_sell1) * df1);
	
	//System.out.println("\n %d curr:%d df0:%f pstate:%f add:%f cash[t]:%f",t,curr,df0,pstate,pstate * bio.getRev()[curr][1][1][0][sea] * df0,mkv.getCash()[t]);
	//System.out.println(" %f ",pid * (- p_loss + p_sell1) * df1);
	//System.out.println("\n curr:%d pid:%f  p_loss:%f  p_sell1:%f  cash[t]:%f",curr,pid,p_loss,p_sell1,mkv.getCash()[t]);
	//hlp += pstate * (bio.getRev()[curr][1][1][0][sea] * df0 + pid * (- p_loss + p_sell1) * df1);
	//System.out.println("\nt:%d mim:%d curr:%d pstate:%f rev:%f %f (%f)",t,mim,curr,pstate,bio.getRev()[curr][1][1][0][sea],pstate*bio.getRev()[curr][1][1][0][sea],hlp);
	
					mkv.getRev()[t] += pstate * bio.getRev()[curr][1][1][0][sea+w] * df0;
					mkv.getMilkyield()[mim][t] += pstate * bio.getMilkyield()[curr][1][1][0][sea];
					mkv.getMilkyield_cnt()[mim][t] += pstate;
					mkv.getFeedcost()[1][t] += pstate * bio.getDmi()[curr][1][1][0][sea] * bio.getP_tmrdmlac() * df0;
					mkv.getDmi()[mim][t] += pstate * bio.getDmi()[curr][lact][mim][0][sea];
					mkv.getDmi_cnt()[mim][t] += pstate;
					mkv.getMilking()[t] += pstate;								//fraction cows milking at beginning of stage
					mkv.getCowsales()[t] += pstate * pid * p_sell1 * df1;		//cows sold at end of stage
					mkv.getLoss()[t] += pstate * pid * p_loss * df1;			//loss due to involuntary culling
					mkv.getOpen()[curr][t] += pstate;							//fraction open at beginning of stage
					mkv.getDim()[t] += pstate * (Constants.STAGE / 2);					//days in milk
					mkv.getLaborcost()[t] += pstate * bio.getTimecows()[curr][1][1][0][sea] / 60 * bio.getP_labor() * df0;
					mkv.getLabortime()[t] += pstate * bio.getTimecows()[curr][1][1][0][sea];
					mkv.getVetcost()[t] += pstate * bio.getVetcost()[1] * df0;		//veterinary cost
					mkv.getCalved()[t] += pstate;								//cows calving
					mkv.getPind()[curr][t] += pstate * pid;
					mkv.getPind_mim()[mim][t] += pstate * pid;
					mkv.getStatus()[mim][0] += pstate;
					mkv.getLact()[1][t] += pstate;								//1st lactation
					mkv.getDysopen()[t] += pstate * (double)(mim);
					mkv.getOthercost()[t] += pstate * bio.getP_other_lact() * Constants.STAGE * df0;
					if(rpo.getCow()[curr][1][1][0][sea+w] != -Constants.X) mkv.getRpo()[t] += pstate * rpo.getCow()[curr][1][1][0][sea+w];
					if(rpo.getCow()[curr][1][1][0][sea+w] != -Constants.X) mkv.getRpo_cnt()[t] += pstate;
					if(rpo.getCow()[curr][1][1][0][sea+w] != -Constants.X) mkv.getPind_rpo()[curr][1][1][0] += pstate * pid * rpo.getCow()[curr][1][1][0][sea+w];
					if(rpo.getCow()[curr][1][1][0][sea+w] != -Constants.X) mkv.getPind_rpo_cnt()[curr][1][1][0] += pstate * pid;
					if(rpo.getCow()[curr][lact][mim][0][sea+w] > 0) mkv.getValue()[t] += pstate * (rpo.getCow()[curr][1][1][0][sea+w] + bio.getSel()[1][1][sea+w]);	//although we don't voluntary sell a heifer
					else mkv.getValue()[t] += pstate * p_sell0;
					mkv.getValue_cnt()[t] += pstate;
					mkv.getHz_open()[1][0] += pstate;							//pind open: at end of month
					mkv.getHz_open()[1][2] += pstate * pid;						//pind open: at end of month
				}
			}
	//*** WHOLE HERD VARIABLES **************************************************************
			if(all == 1)
			{
				for(i = 1; i <= Constants.MIM; i++) mkv.getCwt()[t] += mkv.getMilkyield()[i][t] / 100;
				for(i = 1; i <= Constants.MIM; i++) mkv.getMilksales()[t] += mkv.getMilkyield()[i][t] * bio.getP_milk()[sea+w] * df0;
				mkv.getHeiferpurchasecost()[t] = mkv.getHeifers_purchased()[t] * p_heifer * df0;
				mkv.getCalfsales()[t] = mkv.getCalves_born()[t] * bio.getP_calf() * df0;
				for(curr = 1; curr <= Constants.CURR; curr++) mkv.getPresent()[curr][t] = mkv.getOpen()[curr][t] + mkv.getPreg()[curr][t];
				mkv.getFixedcost_labor()[t] = bio.getP_fixed_labor() * Constants.STAGE * df0;
				mkv.getFixedcost_other()[t] = bio.getP_fixed_other() * Constants.STAGE * df0;
				mkv.getCash()[t] -= (mkv.getFixedcost_labor()[t] + mkv.getFixedcost_other()[t]);
	
				mkv.getTotrevenue()[t] = mkv.getMilksales()[t] + mkv.getCowsales()[t] + mkv.getCalfsales()[t];
				mkv.getTotcost()[t] = mkv.getFeedcost()[0][t] + mkv.getFeedcost()[1][t] + mkv.getBreedingsupplycost()[t] 
								+ mkv.getHeiferpurchasecost()[t] + mkv.getLaborcost()[t] + mkv.getVetcost()[t] + mkv.getLoss()[t] 
								+ mkv.getOthercost()[t] + mkv.getFixedcost_labor()[t] + mkv.getFixedcost_other()[t];
	
				if((int)t % (int)((double)maxstages / 5) == 0 || t == 2 || t == maxstages)
				{
					System.out.println(50+(double)t/(double)maxstages*50);
					Utility.write_msg_file(15,50+(double)t/(double)maxstages*50);
				}
	
			//	System.out.println("\n t:%d %d %d %f",t,(int)maxstages,(int)((double)maxstages / 10),(int)maxstages % (int)((double)maxstages / 10));
			}
	
	
	//		for(k = 1; k <= 2; k++)	//Constants.MIM
	//		{
	//			mkvsum = 0.0;
	//			for(i = 0; i <= Constants.CURR; i++)
	//				for(j = 0; j <= Constants.LACT; j++)
	//					mkvsum += mkv.getCow()[i][j][k][0][ii];	//ending stage
	//System.out.println("\n t:%d mim:%d  %f  %f",t,k,mkvsum,mkvsum+mkv.getheiferpool[t]);
	//		}
	//stop(0);
	
	
			mkvsum = 0.0;
			for(i = 0; i <= Constants.CURR; i++)
				for(j = 0; j <= bio.getMaxlact(); j++)
					for(k = 0; k <= Constants.MIM; k++)
						for(m = 0; m <= Constants.MPREG; m++)
							mkvsum += mkv.getCow()[i][j][k][m][ii];	//entering stage
	//System.out.println("\n%d mkvsumii:%f  at start",t,mkvsum);
			mkvsum = 0.0;
			for(i = 0; i <= Constants.CURR; i++)
				for(j = 0; j <= bio.getMaxlact(); j++)
					for(k = 0; k <= Constants.MIM; k++)
						for(m = 0; m <= Constants.MPREG; m++)
							mkvsum += mkv.getCow()[i][j][k][m][jj];	//ending stage
	//System.out.println("\n%d mkvsumjj:%f  at end",t,mkvsum);
	//stop(0);
	//f(t <= 3) stop(0);
	
			jj = ii;								//swap ii(this stage) and jj(next stage)
			if(ii == 0) ii = 1; else ii = 0;
			sea++; if(sea > Constants.SEA) sea = 1;			//season in next stage
			if(t == maxstages)						//store markov chain at stage = maxstage+1 in [2] for writing
			{
				mkv.setFinish_sea(sea);
				for(i = 0; i <= Constants.CURR; i++)
					for(j = 0; j <= bio.getMaxlact(); j++)
						for(k = 0; k <= Constants.MIM; k++)
							for(m = 0; m <= Constants.MPREG; m++)
							{
								mkv.getCow()[i][j][k][m][2] = mkv.getCow()[i][j][k][m][ii];
								mkv.getFs()[i][j][k][2] = mkv.getFs()[i][j][k][ii];
								mkv.getVw()[i][j][k][2] = mkv.getVw()[i][j][k][ii];
							}
				for(i = 0; i <= Constants.AGE; i++)
					for(m = 0; m <= Constants.MPREG; m++)
						mkv.getHeifer()[i][m][2] = mkv.getHeifer()[i][m][ii];
			}
			mkv = InitFunctions.reset_mkv_record(mkv,jj,maxstages);					//reset Pr(next stage) = 0.0
		}//end for(t)
	//	System.out.println("\n MARKOV_CHAIN done ...");
		return;
	}/*** end MARKOV_CHAIN ***/
	
	
	
	
	
	/******************************************************************************
	MARKOV_CHAIN_STATS: Herd statistics calculated with the Markov Chain.
	******************************************************************************/
	@SuppressWarnings("unused")
	public void markov_chain_stats(int maxstages, BioRecord bio, MkvRecord mkv)
	{
		
		int				t = 0, i = 0, s = 0, y = 0, w = 0, sea = 0, curr = 0, lact = 0, mim = 0, mpreg = 0, tt = 0;
		int				cs0 = 0, cs1 = 0, ct0 = 0, ct1 = 0;
		double			years = 0.0, sum = 0.0, df0 = 0.0, df1 = 0.0,
						pv0 = 0.0, pv10 = 0.0, pv0_cnt = 0.0, pv10_cnt = 0.0;
		double repl_open_dim[] = new double[Constants.SEA+1];
		double repl_open[] = new double[Constants.SEA+1];
		double eligible_preg[] = new double[Constants.SEA+1];
		double eligible_keep[] = new double[Constants.SEA+1];
		double eligible_serv[] = new double[Constants.SEA+1];
		
		
		double daysvwp_clv[][] = new double[Constants.LACT+1][Constants.SEA+1];
		double daysvwp_clv_cnt[][] = new double[Constants.LACT+1][Constants.SEA+1];
		double daysfs_clv[][] = new double[Constants.LACT+1][Constants.SEA+1];
		double daysfs_clv_cnt[][] = new double[Constants.LACT+1][Constants.SEA+1];
		double dayslbm_clv[][] = new double[Constants.LACT+1][Constants.SEA+1];
		double dayslbm_clv_cnt[][] = new double[Constants.LACT+1][Constants.SEA+1];
		double daystoconception_clv[][] = new double[Constants.LACT+1][Constants.SEA+1];
		double daystoconception_clv_cnt[][] = new double[Constants.LACT+1][Constants.SEA+1];
		double calvinginterval_clv[][] = new double[Constants.LACT+1][Constants.SEA+1];
		double calvinginterval_clv_cnt[][] = new double[Constants.LACT+1][Constants.SEA+1];
		
		
		double repl_clv[][] = new double[Constants.SEA+1][Constants.YEAR+1];
		double repl_clv_cnt[][] = new double[Constants.SEA+1][Constants.YEAR+1];
		double pind_clv[][] = new double[Constants.SEA+1][Constants.YEAR+1];
		double pind_clv_cnt[][] = new double[Constants.SEA+1][Constants.YEAR+1];
		double 			replopen_cnt = 0.0, replopen_dim = 0.0, replopen_cur = 0.0, replopen_mlk = 0.0,
						replopen_iof = 0.0, replopen_rev = 0.0, replopen_rpo = 0.0, pvifa0 = 0.0, pvifa1 = 0.0;

		
		HRecord h = new HRecord();
		System.out.println("\n MARKOV_CHAIN_STATS started ...");
		//*** initialize local arrays ******************************************************************
		
		for(sea = 0; sea <= Constants.SEA; sea++)
		{
			repl_open_dim[sea] = 0.0;
			repl_open[sea] = 0.0;
			eligible_preg[sea] = 0.0;
			eligible_keep[sea] = 0.0;
			eligible_serv[sea] = 0.0;
	//&
			for(lact = 0; lact <= Constants.LACT; lact++) daysvwp_clv[lact][sea] = 0.0;
			for(lact = 0; lact <= Constants.LACT; lact++) daysvwp_clv_cnt[lact][sea] = 0.0;
			for(lact = 0; lact <= Constants.LACT; lact++) daysfs_clv[lact][sea] = 0.0;
			for(lact = 0; lact <= Constants.LACT; lact++) daysfs_clv_cnt[lact][sea] = 0.0;
			for(lact = 0; lact <= Constants.LACT; lact++) daystoconception_clv[lact][sea] = 0.0;
			for(lact = 0; lact <= Constants.LACT; lact++) daystoconception_clv_cnt[lact][sea] = 0.0;
			for(lact = 0; lact <= Constants.LACT; lact++) dayslbm_clv[lact][sea] = 0.0;
			for(lact = 0; lact <= Constants.LACT; lact++) dayslbm_clv_cnt[lact][sea] = 0.0;
			for(lact = 0; lact <= Constants.LACT; lact++) calvinginterval_clv[lact][sea] = 0.0;
			for(lact = 0; lact <= Constants.LACT; lact++) calvinginterval_clv_cnt[lact][sea] = 0.0;
	
			for(y = 0; y <= Constants.YEAR; y++)
			{
				repl_clv[sea][y] = 0.0;
				repl_clv_cnt[sea][y] = 0.0;
				pind_clv[sea][y] = 0.0;
				pind_clv_cnt[sea][y] = 0.0;
			}
		}
		for(mim = 0; mim <= Constants.MIM; mim++)
		{
			h.getAtrisk()[mim] = 0.0;
			h.getCulled()[mim] = 0.0;
			h.getReplmlkcurr()[mim] = 0.0;
			h.getReplmlkcurr_cnt()[mim] = 0.0;
		}
		for(t = 0; t <= Constants.MAXM; t++)
		{
			h.getClvinterval_cnt()[t] = 0.0;
			for(curr = 0; curr <= Constants.CURR; curr++) h.getDlb_cnt()[curr][t] = 0.0;
			for(curr = 0; curr <= Constants.CURR; curr++) for(lact = 0; lact <= Constants.LACT; lact++) h.getDlbm()[curr][lact][t] = 0.0;
			for(curr = 0; curr <= Constants.CURR; curr++) for(lact = 0; lact <= Constants.LACT; lact++) h.getDlbm_cnt()[curr][lact][t] = 0.0;
			for(mim = 0; mim <= Constants.MIM; mim++)
			{
				h.getInsval()[mim][t] = 0.0;
				h.getInsval_cnt()[mim][t] = 0.0;
				for(curr = 0; curr <= Constants.CURR; curr++)
				{
					h.getReplmlk()[curr][mim][t] = 0.0;
					h.getReplmlk_cnt()[curr][mim][t] = 0.0;
	//&				
					for(lact = 0; lact <= Constants.LACT; lact++) h.getDvwp_cnt()[curr][lact][mim][t] = 0.0;
					for(lact = 0; lact <= Constants.LACT; lact++) h.getDfs_cnt()[curr][lact][mim][t] = 0.0;
				}
			}
		}
	
	
	//*** collect data per year ***************************************************************
		mkv.getCash()[0] = 0.0;	//may have been filled in markov_chain() to monitor total cash sum
		y = 0;
		sea = bio.getStart_sea();
		for (t = 1; t <= maxstages; t++)
		{
			tt++;
		//	if(bio.getAnalysis()== 4) mkv.getherdlife += mkv.getRepl()[t] * (double)(t-1) + mkv.getPind()[t] * (double)t;
		//	not correct because some cows that abort are involuntary culled as well
			if(t <= 12) w = 0; else w = 12;
	
			df0 = Utility.discount(bio.getInterest()/12,(double)tt-1);
			df1 = Utility.discount(bio.getInterest()/12,(double)tt);			//at end of stage: cowsales + loss.
			if(bio.getAnalysis() == 1) df0 = df1 = 1.0;					//steady-state: do not Utility.discount: if discounting is used,
																	//  and there is seasonality, then the starting month has 
			pvifa0 += df0;											//  an effect. With steady-state conditions, the starting
			pvifa1 += df1;											//  month should not matter.
			
			mkv.getHeiferpurchasecost()[0] += bio.getP_heifer()[sea+w] * df0 * mkv.getHeiferpool()[t];
			mkv.getAbortions()[0] += mkv.getAbortions()[t];
			mkv.getAbortions_cnt()[0] += mkv.getAbortions_cnt()[t];
			mkv.getBreedingsupplycost()[0] += mkv.getBreedingsupplycost()[t] * df0;	//drugs + semen
			mkv.getBreedinglaborcost()[0] += mkv.getBreedinglaborcost()[t] * df0;	//labor cost
			mkv.getCalfsales()[0] += mkv.getCalfsales()[t] * df0;
			mkv.getCalved()[0] += mkv.getCalved()[t];
			mkv.getCash()[0] += mkv.getCash()[t];
			mkv.setCash_maxm(mkv.getCash_maxm()+mkv.getCash()[t]);
	
			if(bio.getAnalysis() == 1 && bio.getInterest() <= 0.0001 && t == 12 && Math.abs(mkv.getCash_maxm() - mkv.getCash12()) > 0.01)
		//	if(t == 12)
			{
			//!!markov chain cash/year in steady state may be different from dp cash/year because of bio.getMaxlact() in dp_policy()
			//	System.out.println("\n bio.getAnalysis():%d  interest:%f  fabs:%f",bio.getAnalysis(),bio.getInterest(), fabs(mkv.getCash_maxm() - cash12));
				System.out.println("\n MARKOV_CHAIN_STATS:\n      difference:"+(mkv.getCash_maxm()-mkv.getCash12())+"\n  mkv.getCash_maxm():"+mkv.getCash_maxm()+"\n        DPcash12:"+mkv.getCash12());
				//if(t == 12) stop(0);
				//note: DP_POLICY allows cows to be in lactations > 12th, but MARKOV_CHAIN() does not
				//      this was done to not cull cows necessarily at the end of the 12th lactation
				//      cash_maxm and cas12 should be the same when no cow would want to be in the 12th
				//      lactation, for example when milk yield in the 12th lactation is set to 0%.
			}
	
			mkv.getCowsales()[0] += mkv.getCowsales()[t] * df1;
			mkv.getCwt()[0] += mkv.getCwt()[t];
			mkv.getValue()[0] += mkv.getValue()[t];
			mkv.getValue_cnt()[0] += mkv.getValue_cnt()[t];
	
	
	
	//&		//calvinginterval
			for(lact = 1; lact <= Constants.LACT; lact++)
				for(mim = 1; mim <= Constants.MIM; mim++)
				{
					mkv.getClvinterval()[0][0][t] += mkv.getClvinterval()[lact][mim][t] * (double)mim;
					mkv.getClvinterval()[0][0][0] += mkv.getClvinterval()[lact][mim][t] * (double)mim;
					h.getClvinterval_cnt()[t] +=  mkv.getClvinterval()[lact][mim][t];
					h.getClvinterval_cnt()[0] +=  mkv.getClvinterval()[lact][mim][t];
				}
			if(h.getClvinterval_cnt()[t] > 0.0) mkv.getClvinterval()[0][0][t] = mkv.getClvinterval()[0][0][t] / h.getClvinterval_cnt()[t];
	//&		//new pregnancies, also used for daystoconc
			for(curr = 1; curr <= Constants.CURR; curr++) 
			{
				for(lact = 1; lact <= Constants.LACT; lact++)
				{
					mkv.getNewpreg()[curr][0][0] += mkv.getNewpreg()[curr][lact][t];
					mkv.getNewpreg()[0][lact][0] += mkv.getNewpreg()[curr][lact][t];
					mkv.getNewpreg()[0][0][t] += mkv.getNewpreg()[curr][lact][t];
					mkv.getNewpreg()[0][0][0] += mkv.getNewpreg()[curr][lact][t];
				}
			}
	//&
			for(curr = 1; curr <= Constants.CURR; curr++)
			{
				for(mim = 1; mim <= Constants.MIM; mim++)
				{
					for(lact = 1; lact <= Constants.LACT; lact++)
					{
						mkv.getDaystoconc()[curr][0][0][0] += mkv.getDaystoconc()[curr][lact][mim][t] * (double)mim * Constants.STAGE;//numerator
						mkv.getDaystoconc()[0][0][0][t] += mkv.getDaystoconc()[curr][lact][mim][t] * (double)mim * Constants.STAGE;	//numerator
						mkv.getDaystoconc()[0][0][0][0] += mkv.getDaystoconc()[curr][lact][mim][t] * (double)mim * Constants.STAGE;	//numerator
	
						mkv.getDvwp()[curr][0][0][0] += mkv.getDvwp()[curr][lact][mim][t] * (double)mim * Constants.STAGE;
						mkv.getDvwp()[0][0][0][t] += mkv.getDvwp()[curr][lact][mim][t] * (double)mim * Constants.STAGE;
						mkv.getDvwp()[0][0][0][0] += mkv.getDvwp()[curr][lact][mim][t] * (double)mim * Constants.STAGE;
						h.getDvwp_cnt()[curr][0][0][0] += mkv.getDvwp()[curr][lact][mim][t];
						h.getDvwp_cnt()[0][0][0][t] += mkv.getDvwp()[curr][lact][mim][t];
						h.getDvwp_cnt()[0][0][0][0] += mkv.getDvwp()[curr][lact][mim][t];
	
						mkv.getDfs()[curr][0][0][0] += mkv.getDfs()[curr][lact][mim][t] * (double)mim * Constants.STAGE;
						mkv.getDfs()[0][0][0][t] += mkv.getDfs()[curr][lact][mim][t] * (double)mim * Constants.STAGE;
						mkv.getDfs()[0][0][0][0] += mkv.getDfs()[curr][lact][mim][t] * (double)mim * Constants.STAGE;
						h.getDfs_cnt()[curr][0][0][0] += mkv.getDfs()[curr][lact][mim][t];
						h.getDfs_cnt()[0][0][0][t] += mkv.getDfs()[curr][lact][mim][t];
						h.getDfs_cnt()[0][0][0][0] += mkv.getDfs()[curr][lact][mim][t];
	
					//	System.out.println("\n t:%d curr:%d mim:%d %f %f",t,curr,mim,mkv.getDaystoconc()[curr][0][0],mkv.getDaystoconc()[curr][0][0]/mkv.getNewpreg()[curr][0]);
					}
					mkv.getDysopencost()[curr][0][0] += mkv.getDysopencost()[curr][mim][t];
					mkv.getDysopencost()[0][mim][0] += mkv.getDysopencost()[curr][mim][t];
					mkv.getDysopencost()[0][0][0] += mkv.getDysopencost()[curr][mim][t];
	
					mkv.getDysopencost_cnt()[curr][0][0] += mkv.getDysopencost_cnt()[curr][mim][t];
					mkv.getDysopencost_cnt()[0][mim][0] += mkv.getDysopencost_cnt()[curr][mim][t];
					mkv.getDysopencost_cnt()[0][0][0] += mkv.getDysopencost_cnt()[curr][mim][t];
	
				//	if(mkv.getDysopencost_cnt()[0][mim][0] > 0)
				//	System.out.println("\n t:%d curr:%d mim:%d %f %f",t,curr,mim,mkv.getDysopencost()[0][mim][0],mkv.getDysopencost()[0][mim][0]/mkv.getDysopencost_cnt()[0][mim][0]);
	
					mkv.getDyspregcost()[curr][0][0] += mkv.getDyspregcost()[curr][mim][t];
					mkv.getDyspregcost()[0][mim][0] += mkv.getDyspregcost()[curr][mim][t];
					mkv.getDyspregcost()[0][0][0] += mkv.getDyspregcost()[curr][mim][t];
	
					mkv.getDyspregcost_cnt()[curr][0][0] += mkv.getDyspregcost_cnt()[curr][mim][t];
					mkv.getDyspregcost_cnt()[0][mim][0] += mkv.getDyspregcost_cnt()[curr][mim][t];
					mkv.getDyspregcost_cnt()[0][0][0] += mkv.getDyspregcost_cnt()[curr][mim][t];
				}
				//	denominator is newpreg[curr][0] or newpreg[0][0]
			}
	//&
			if(mkv.getNewpreg()[0][0][t] > 0.0) mkv.getDaystoconc()[0][0][0][t] = mkv.getDaystoconc()[0][0][0][t] / mkv.getNewpreg()[0][0][t];
			if(h.getDvwp_cnt()[0][0][0][t] > 0.0) mkv.getDvwp()[0][0][0][t] = mkv.getDvwp()[0][0][0][t] / h.getDvwp_cnt()[0][0][0][t];
			if(h.getDfs_cnt()[0][0][0][t] > 0.0) mkv.getDfs()[0][0][0][t] = mkv.getDfs()[0][0][0][t] / h.getDfs_cnt()[0][0][0][t];
	
	//&
			for(lact = 1; lact <= Constants.LACT; lact++)
			{
				for(curr = 1; curr <= Constants.CURR; curr++)
				{
					mim = Constants.MIM + 1;
					do
					{
						mim = mim - 1;
						h.getDlbm()[curr][lact][t] = mkv.getDlbm()[curr][lact][mim][t] * (double)mim * Constants.STAGE;
						h.getDlbm_cnt()[curr][lact][t] = mkv.getDlbm()[curr][lact][mim][t];
					
					//	System.out.println("\nt:%d lact:%d curr:%d mim:%d  mkv.getDlbm():%f  h.getDlbm():%f",
					//		t,lact,curr,mim,mkv.getDlbm()[curr][lact][mim][t],h.getDlbm()[curr][lact][t]);
					//	if(h.getDlbm()[curr][lact][t] >= 740 || mim > Constants.MIM) stop(0);
					}
					while(mkv.getDlbm()[curr][lact][mim][t] < 0.000000001 && mim > 1);
				}
			}
			for(lact = 1; lact <= Constants.LACT; lact++)	
				for(curr = 1; curr <= Constants.CURR; curr++)
				{
					h.getDlbm()[0][0][t] += h.getDlbm()[curr][lact][t];
					h.getDlbm()[0][0][0] += h.getDlbm()[curr][lact][t];
					h.getDlbm_cnt()[0][0][t] += h.getDlbm_cnt()[curr][lact][t];
					h.getDlbm_cnt()[0][0][0] += h.getDlbm_cnt()[curr][lact][t];
				}
			if(h.getDlbm_cnt()[0][0][t] > 0.0) mkv.getDlbm()[0][0][0][t] = h.getDlbm()[0][0][t] / h.getDlbm_cnt()[0][0][t];
	
	
	
			mkv.getDysopen()[0] += mkv.getDysopen()[t];
			mkv.getFeedcost()[0][0] += mkv.getFeedcost()[0][t] * df0;			//dry cows
			mkv.getFeedcost()[1][0] += mkv.getFeedcost()[1][t] * df0;			//lactating cows
			mkv.getFixedcost_labor()[0] += mkv.getFixedcost_labor()[t] * df0;
			mkv.getFixedcost_other()[0] += mkv.getFixedcost_other()[t] * df0;
			mkv.getMilking()[0] += mkv.getMilking()[t];
			for(mim = 1; mim <= Constants.MIM; mim++) mkv.getMilkyield()[0][t] += mkv.getMilkyield()[mim][t];
			mkv.getMilkyield()[0][0] += mkv.getMilkyield()[0][t];
			mkv.getMilksales()[0] += mkv.getMilksales()[t] * df0;
			mkv.getLaborcost()[0] += mkv.getLaborcost()[t] * df0;
			mkv.getLoss()[0] += mkv.getLoss()[t] * df1;
			mkv.getOthercost()[0] += mkv.getOthercost()[t] * df0;
	
			for(curr = 1; curr <= Constants.CURR; curr++) 
			{
				mkv.getPreg()[0][0] += mkv.getPreg()[curr][t];
				mkv.getPreg()[curr][0] += mkv.getPreg()[curr][t];
				mkv.getPreg()[0][t] += mkv.getPreg()[curr][t];
	
				mkv.getKeep()[0][0] += mkv.getKeep()[curr][t];
				mkv.getKeep()[curr][0] += mkv.getKeep()[curr][t];
				mkv.getKeep()[0][t] += mkv.getKeep()[curr][t];
	
				mkv.getOpen()[0][0] += mkv.getOpen()[curr][t];
				mkv.getOpen()[curr][0] += mkv.getOpen()[curr][t];
				mkv.getOpen()[0][t] += mkv.getOpen()[curr][t];
	
				mkv.getPind()[0][0] += mkv.getPind()[curr][t];
				mkv.getPind()[curr][0] += mkv.getPind()[curr][t];
				mkv.getPind()[0][t] += mkv.getPind()[curr][t];
			
				mkv.getPresent()[0][0] += mkv.getPresent()[curr][t];
				mkv.getPresent()[curr][0] += mkv.getPresent()[curr][t];
				mkv.getPresent()[0][t] += mkv.getPresent()[curr][t];
			}
			if(mkv.getPresent()[0][t] > 0.0) for(curr = 1; curr <= Constants.CURR; curr++)
				mkv.getSmy()[t] += mkv.getPresent()[curr][t] / mkv.getPresent()[0][t] * bio.getSmy()[curr];	//average smy per t
	
			for(curr = 1; curr <= Constants.CURR; curr++) mkv.getRepl()[0][0] += mkv.getRepl()[curr][t];
			for(curr = 1; curr <= Constants.CURR; curr++) mkv.getRepl()[curr][0] += mkv.getRepl()[curr][t];
			for(curr = 1; curr <= Constants.CURR; curr++) mkv.getRepl()[0][t] += mkv.getRepl()[curr][t];
	
			mkv.getRpo()[0] += mkv.getRpo()[t];
			mkv.getRpo_cnt()[0] += mkv.getRpo_cnt()[t];
			mkv.getVetcost()[0] += mkv.getVetcost()[t] * df0;
	
	//if(t == 14)
	//{
	//	for(curr = 1; curr <= Constants.CURR; curr++)
	//	{
	//		System.out.println("\n %2d",curr);
	//		for(i = 1; i <= t; i++)
	//		{
	//			System.out.println(" %.5f",mkv.getPresent()[curr][i]);
	//		//	System.out.println(" %.5f",mkv.getPreg()[curr][i]);
	//		}
	//	}
	//}
	//if(t == 24) stop(0);
	
	
	//observation:
	//if delay is an option, all high cows are delayed and low cows are bred -.get low cows calve again and then because
	//of regression towards the mean, start producing a bit more. the total is that above 100 yield is reached.
	//if delay is temporily, all high are delayed a while, so the low producers calve first (yeidl goes up) but 
	//shortly after, all high producers start calving -.get yield goes down because of regression towards the mean
	//-.get is should all cancel out.
	//at t=12, mkv_stats shows a higher milk yield (low producers calve + regression towards the mean, sounds ok)
	//however, dp_policy does not show that, why not?  or are sea and t inverse?
	
	//		System.out.println("\n t:%d df0:%f  milking:%f yield:%.3f cash:%f   %f",t,df0,mkv.getMilking()[t],mkv.getMilkyield()[t],mkv.getCash()[t],mkv.getMilksales()[0]  + mkv.getCalfsales()[0]
	//				- mkv.getFeedcost()[0][0] - mkv.getFeedcost()[1][0] - mkv.getBreedingsupplycost()[0] 
	//				- mkv.getHeiferpurchasecost()[0] - mkv.getVetcost()[0] - mkv.getLaborcost()[0] - mkv.getOthercost()[0] 
	//				- mkv.getFixedcost_labor()[0] - mkv.getFixedcost_other()[0]
	//				+ mkv.getCowsales()[0] - mkv.getLoss()[0]);
	//
	//		System.out.println(" milksales:%f",mkv.getMilksales()[0]); 
	
			for(mim = 1; mim <= Constants.MIM; mim++)	
			{
				eligible_keep[0] += mkv.getEligible_keep()[mim][t];
				eligible_preg[0] += mkv.getEligible_preg()[mim][t];
				eligible_serv[0] += mkv.getEligible_serv()[mim][t];
	
				//seasonal variables
				eligible_keep[sea] += mkv.getEligible_keep()[mim][t];
				eligible_preg[sea] += mkv.getEligible_preg()[mim][t];
				eligible_serv[sea] += mkv.getEligible_serv()[mim][t];
			}
			for(mim = 1; mim <= Constants.MIM; mim++)	
			{
				mkv.getEligible_bred()[0][t] += mkv.getEligible_bred()[mim][t];
				mkv.getEligible_keep()[0][t] += mkv.getEligible_keep()[mim][t];
				mkv.getEligible_preg()[0][t] += mkv.getEligible_preg()[mim][t];
				mkv.getEligible_serv()[0][t] += mkv.getEligible_serv()[mim][t];
			}
	
	
			for(mim = 1; mim <= Constants.MIM; mim++)	for(curr = 1; curr <= Constants.CURR; curr++)
			{
				//probability
				mkv.getReplopen_cnt()[0][t][0]   += mkv.getRepl_open_cnt()[curr][mim][t];
				mkv.getReplopen_cnt()[mim][0][0] += mkv.getRepl_open_cnt()[curr][mim][t];
				mkv.getReplopen_cnt()[0][0][0]   += mkv.getRepl_open_cnt()[curr][mim][t];
				replopen_cnt                 += mkv.getRepl_open_cnt()[curr][mim][t];
	
				//dim at voluntary culling open cows
				mkv.getReplopen_dim()[0][t][0]   += (double)(mim - 1) * mkv.getRepl_open_cnt()[curr][mim][t] * Constants.STAGE;
				mkv.getReplopen_dim()[mim][0][0] += (double)(mim - 1) * mkv.getRepl_open_cnt()[curr][mim][t] * Constants.STAGE;
				mkv.getReplopen_dim()[0][0][0]   += (double)(mim - 1) * mkv.getRepl_open_cnt()[curr][mim][t] * Constants.STAGE;
				replopen_dim                 += (double)(mim - 1) * mkv.getRepl_open_cnt()[curr][mim][t] * Constants.STAGE;
	
				//curr at voluntary culling open cows
				mkv.getReplopen_cur()[0][t][0]   += (double)curr * mkv.getRepl_open_cnt()[curr][mim][t];
				mkv.getReplopen_cur()[mim][0][0] += (double)curr * mkv.getRepl_open_cnt()[curr][mim][t];
				mkv.getReplopen_cur()[0][0][0]   += (double)curr * mkv.getRepl_open_cnt()[curr][mim][t];
				replopen_cur                 += (double)curr * mkv.getRepl_open_cnt()[curr][mim][t];
	
				//milk yield at voluntary culling open cows
				mkv.getReplopen_mlk()[0][t][0]   += mkv.getRepl_open_mlk()[curr][mim][t] / Constants.STAGE;
				mkv.getReplopen_mlk()[mim][0][0] += mkv.getRepl_open_mlk()[curr][mim][t] / Constants.STAGE;
				mkv.getReplopen_mlk()[0][0][0]   += mkv.getRepl_open_mlk()[curr][mim][t] / Constants.STAGE;
				replopen_mlk                 += mkv.getRepl_open_mlk()[curr][mim][t] / Constants.STAGE;
	
				//rpo at voluntary culling open cows
				mkv.getReplopen_rpo()[0][t][0]   += mkv.getRepl_open_rpo()[curr][mim][t];
				mkv.getReplopen_rpo()[mim][0][0] += mkv.getRepl_open_rpo()[curr][mim][t];
				mkv.getReplopen_rpo()[0][0][0]   += mkv.getRepl_open_rpo()[curr][mim][t];
				replopen_rpo                 += mkv.getRepl_open_rpo()[curr][mim][t];
	
				//iofc at voluntary culling open cows
				mkv.getReplopen_iof()[0][t][0]   += mkv.getRepl_open_iof()[curr][mim][t] / Constants.STAGE;
				mkv.getReplopen_iof()[mim][0][0] += mkv.getRepl_open_iof()[curr][mim][t] / Constants.STAGE;
				mkv.getReplopen_iof()[0][0][0]   += mkv.getRepl_open_iof()[curr][mim][t] / Constants.STAGE;
				replopen_iof                 += mkv.getRepl_open_iof()[curr][mim][t] / Constants.STAGE;
	
				//margin at voluntary culling open cows
				mkv.getReplopen_rev()[0][t][0]   += mkv.getRepl_open_rev()[curr][mim][t] / Constants.STAGE;
				mkv.getReplopen_rev()[mim][0][0] += mkv.getRepl_open_rev()[curr][mim][t] / Constants.STAGE;
				mkv.getReplopen_rev()[0][0][0]   += mkv.getRepl_open_rev()[curr][mim][t] / Constants.STAGE;
				replopen_rev                 += mkv.getRepl_open_rev()[curr][mim][t] / Constants.STAGE;
	
			//	System.out.println("\nt:%d mim:%d curr:%d",t,mim,curr);
			//	System.out.println("\ndim  %.15f  %.15f %f",mkv.getReplopen_dim()[0][0][0],mkv.getReplopen_cnt()[0][0][0],mkv.getReplopen_dim()[0][0][0]/mkv.getReplopen_cnt()[0][0][0]);
			//	System.out.println("\ncur  %.15f  %.15f %f",mkv.getReplopen_cur()[0][0][0],mkv.getReplopen_cnt()[0][0][0],mkv.getReplopen_cur()[0][0][0]/mkv.getReplopen_cnt()[0][0][0]);
			//	System.out.println("\nmlk  %.15f  %.15f %f",mkv.getReplopen_mlk()[0][0][0],mkv.getReplopen_cnt()[0][0][0],mkv.getReplopen_mlk()[0][0][0]/mkv.getReplopen_cnt()[0][0][0]);
			//	System.out.println("\nrev  %.15f  %.15f %f",mkv.getReplopen_rev()[0][0][0],mkv.getReplopen_cnt()[0][0][0],mkv.getReplopen_rev()[0][0][0]/mkv.getReplopen_cnt()[0][0][0]);
			//	System.out.println("\niof  %.15f  %.15f %f",mkv.getReplopen_iof()[0][0][0],mkv.getReplopen_cnt()[0][0][0],mkv.getReplopen_iof()[0][0][0]/mkv.getReplopen_cnt()[0][0][0]);
			//	System.out.println("\nrpo  %.15f  %.15f %f",mkv.getReplopen_rpo()[0][0][0],mkv.getReplopen_cnt()[0][0][0],mkv.getReplopen_rpo()[0][0][0]/mkv.getReplopen_cnt()[0][0][0]);
			//	stop(0);
			}
	
			//pregnancyvalue
			for(curr = 1; curr <= Constants.CURR; curr++) for(lact = 1; lact <= bio.getMaxlact(); lact++) for(mim = 1; mim <= Constants.MIM; mim++) for(mpreg = 0; mpreg <= Constants.MPREG; mpreg++)	
			{
				mkv.getPv_curr()[curr][mpreg] += mkv.getPregval()[curr][lact][mim][mpreg][t];
				mkv.getPv_curr_cnt()[curr][mpreg] += mkv.getPregval_cnt()[curr][lact][mim][mpreg][t];
				if(mpreg >= 1) mkv.getPv_curr()[curr][10] += mkv.getPregval()[curr][lact][mim][mpreg][t];
				if(mpreg >= 1) mkv.getPv_curr_cnt()[curr][10] += mkv.getPregval_cnt()[curr][lact][mim][mpreg][t];
	
				mkv.getPv_lact()[lact][mpreg] += mkv.getPregval()[curr][lact][mim][mpreg][t];
				mkv.getPv_lact_cnt()[lact][mpreg] += mkv.getPregval_cnt()[curr][lact][mim][mpreg][t];
				if(mpreg >= 1) mkv.getPv_lact()[lact][10] += mkv.getPregval()[curr][lact][mim][mpreg][t];
				if(mpreg >= 1) mkv.getPv_lact_cnt()[lact][10] += mkv.getPregval_cnt()[curr][lact][mim][mpreg][t];
	
				mkv.getPv_mim()[mim][mpreg] += mkv.getPregval()[curr][lact][mim][mpreg][t];
				mkv.getPv_mim_cnt()[mim][mpreg] += mkv.getPregval_cnt()[curr][lact][mim][mpreg][t];
				if(mpreg >= 1) mkv.getPv_mim()[mim][10] += mkv.getPregval()[curr][lact][mim][mpreg][t];
				if(mpreg >= 1) mkv.getPv_mim_cnt()[mim][10] += mkv.getPregval_cnt()[curr][lact][mim][mpreg][t];
	
				mkv.getPv_mpreg()[mpreg] += mkv.getPregval()[curr][lact][mim][mpreg][t];
				mkv.getPv_mpreg_cnt()[mpreg] += mkv.getPregval_cnt()[curr][lact][mim][mpreg][t];
				if(mpreg >= 1) mkv.getPv_mpreg()[10] += mkv.getPregval()[curr][lact][mim][mpreg][t];
				if(mpreg >= 1) mkv.getPv_mpreg_cnt()[10] += mkv.getPregval_cnt()[curr][lact][mim][mpreg][t];
	
				mkv.getPv_t()[t][mpreg] += mkv.getPregval()[curr][lact][mim][mpreg][t];
				mkv.getPv_t_cnt()[t][mpreg] += mkv.getPregval_cnt()[curr][lact][mim][mpreg][t];
				if(mpreg >= 1) mkv.getPv_t()[t][10] += mkv.getPregval()[curr][lact][mim][mpreg][t];
				if(mpreg >= 1) mkv.getPv_t_cnt()[t][10] += mkv.getPregval_cnt()[curr][lact][mim][mpreg][t];
	
				mkv.getPv_curr_lact()[curr][lact][mpreg] += mkv.getPregval()[curr][lact][mim][mpreg][t];
				mkv.getPv_curr_lact_cnt()[curr][lact][mpreg] += mkv.getPregval_cnt()[curr][lact][mim][mpreg][t];
				if(mpreg >= 1) mkv.getPv_curr_lact()[curr][lact][10] += mkv.getPregval()[curr][lact][mim][mpreg][t];
				if(mpreg >= 1) mkv.getPv_curr_lact_cnt()[curr][lact][10] += mkv.getPregval_cnt()[curr][lact][mim][mpreg][t];
	
				if(mpreg >= 1) mkv.getPv_mpreg()[10] += mkv.getPregval()[curr][lact][mim][mpreg][t];
				if(mpreg >= 1) mkv.getPv_mpreg_cnt()[10] += mkv.getPregval_cnt()[curr][lact][mim][mpreg][t];
	
				if(mpreg == 0) pv0 += mkv.getPregval()[curr][lact][mim][mpreg][t];
				if(mpreg == 0) pv0_cnt += mkv.getPregval_cnt()[curr][lact][mim][mpreg][t];
				if(mpreg >= 1) pv10 += mkv.getPregval()[curr][lact][mim][mpreg][t];
				if(mpreg >= 1) pv10_cnt += mkv.getPregval_cnt()[curr][lact][mim][mpreg][t];
			}
	
			//insemination value
			for(mim = 1; mim <= Constants.MIM; mim++) 
			{
				h.getInsval()[mim][0] += mkv.getInsval()[mim][t];
				h.getInsval()[0][t] += mkv.getInsval()[mim][t];
				h.getInsval()[0][0] += mkv.getInsval()[mim][t];
				h.getInsval_cnt()[mim][0] += mkv.getEligible_keep()[mim][t] + mkv.getEligible_serv()[mim][t];
				h.getInsval_cnt()[0][t] += mkv.getEligible_keep()[mim][t] + mkv.getEligible_serv()[mim][t];
				h.getInsval_cnt()[0][0] += mkv.getEligible_keep()[mim][t] + mkv.getEligible_serv()[mim][t];
			}
	
	
			//cullprice
			if(mkv.getPind()[0][t]+mkv.getRepl()[0][t] > 0.0)
				mkv.getCullprice()[t] = mkv.getCowsales()[t]/(mkv.getPind()[0][t]+mkv.getRepl()[0][t]);
			//days in milk
			mkv.getDim()[0] += mkv.getDim()[t];
			if(mkv.getMilking()[t] > 0) mkv.getDim_()[t] = mkv.getDim()[t] / mkv.getMilking()[t];
	
	
	
	//*** calculate statistics per year *******************************************************
			if(t == 12 || t == 24 || t == 36 || t == maxstages)	
			{
				if(t == 12) y = 1; else if(t == 24) y = 2; else if(t == 36) y = 3; else y = 4;
				if(t == 12 || t == 24 || t == 36) years = 1.0; else years = ((double)maxstages - 36)/12;
				tt = 0;	//time counter for discounting of annual statistics: EAA per year
				//add statistics per stall != per cow, when delay occurs
	
			//	if(bio.getInterest() > 0)
			//		mkv.getannuity_cyr[0] = Constants.SEA * mkv.getCash()[0] * bio.getInterest()/Constants.SEA / (1 - pow(1 + bio.getInterest()/Constants.SEA,-Constants.MAXM));
				//Rogers[702]3454 and Battles+Thompson[Fundamentals of Agribusiness Finance]27
				//cash flows are monthly discounted -.get the annuity is first calculated per month.
				//because the annuity is equal in every month, the annuity per year is 12 * monthly annuity
				//note: this is the annual annuity per stall, not per cow
	
	//System.out.println("\n %.2f",mkv.getMilksales()[0]);
	//System.out.println("\n %.2f",mkv.getCalfsales()[0]);
	//System.out.println("\n %.2f",mkv.getFeedcost()[0][0] + mkv.getFeedcost()[1][0]);
	//System.out.println("\n %.2f",mkv.getBreedingsupplycost()[0]);
	//System.out.println("\n %.2f",mkv.getHeiferpurchasecost()[0]);
	//System.out.println("\n %.2f",mkv.getVetcost()[0]);
	//System.out.println("\n %.2f",mkv.getLaborcost()[0]);
	//System.out.println("\n %.2f",mkv.getOthercost()[0]);
	//System.out.println("\n %.2f",mkv.getFixedcost_labor()[0]);
	//System.out.println("\n %.2f",mkv.getFixedcost_other()[0]);
	//System.out.println("\n %.2f",mkv.getCowsales()[0]);
	//System.out.println("\n %.2f",mkv.getLoss()[0]);
	//stop(0);
	
	
				mkv.getMilksales()[0] = mkv.getMilksales()[0] / pvifa0 * 12;
				mkv.getCalfsales()[0] = mkv.getCalfsales()[0] / pvifa0 * 12;
				mkv.getFeedcost()[0][0] = mkv.getFeedcost()[0][0] / pvifa0 * 12;
				mkv.getFeedcost()[1][0] = mkv.getFeedcost()[1][0] / pvifa0 * 12;
				mkv.getBreedingsupplycost()[0] = mkv.getBreedingsupplycost()[0] / pvifa0 * 12;
				mkv.getHeiferpurchasecost()[0] = mkv.getHeiferpurchasecost()[0] / pvifa0 * 12;
				mkv.getVetcost()[0] = mkv.getVetcost()[0] / pvifa0 * 12;
				mkv.getLaborcost()[0] = mkv.getLaborcost()[0] / pvifa0 * 12;
				mkv.getOthercost()[0] = mkv.getOthercost()[0] / pvifa0 * 12;
				mkv.getFixedcost_labor()[0] = mkv.getFixedcost_labor()[0] / pvifa0 * 12;
				mkv.getFixedcost_other()[0] = mkv.getFixedcost_other()[0] / pvifa0 * 12;
				mkv.getCowsales()[0] = mkv.getCowsales()[0] / pvifa1 * 12;
				mkv.getLoss()[0] = mkv.getLoss()[0] / pvifa1 * 12;
	
	//System.out.println("\n %.2f",mkv.getMilksales()[0]);
	//System.out.println("\n %.2f",mkv.getCalfsales()[0]);
	//System.out.println("\n %.2f",mkv.getFeedcost()[0][0] + mkv.getFeedcost()[1][0]);
	//System.out.println("\n %.2f",mkv.getBreedingsupplycost()[0]);
	//System.out.println("\n %.2f",mkv.getHeiferpurchasecost()[0]);
	//System.out.println("\n %.2f",mkv.getVetcost()[0]);
	//System.out.println("\n %.2f",mkv.getLaborcost()[0]);
	//System.out.println("\n %.2f",mkv.getOthercost()[0]);
	//System.out.println("\n %.2f",mkv.getFixedcost_labor()[0]);
	//System.out.println("\n %.2f",mkv.getFixedcost_other()[0]);
	//System.out.println("\n %.2f",mkv.getCowsales()[0]);
	//System.out.println("\n %.2f",mkv.getLoss()[0]);
	//stop(0);
	
				mkv.getEaa_cyr()[y] = (mkv.getMilksales()[0]  + mkv.getCalfsales()[0]
					- mkv.getFeedcost()[0][0] - mkv.getFeedcost()[1][0] - mkv.getBreedingsupplycost()[0] 
					- mkv.getHeiferpurchasecost()[0] - mkv.getVetcost()[0] - mkv.getLaborcost()[0] - mkv.getOthercost()[0] 
					- mkv.getFixedcost_labor()[0] - mkv.getFixedcost_other()[0]
					+ mkv.getCowsales()[0] - mkv.getLoss()[0]) / years;
				
				System.out.println("MKV CWT : "+mkv.getCwt()[0]);
				mkv.setCwt_final(mkv.getCwt()[0]);
	
				if(mkv.getCwt()[0] > 0.0001)
				mkv.getEaa_cwt()[y] = (mkv.getMilksales()[0] + mkv.getCalfsales()[0]
					- mkv.getFeedcost()[0][0] - mkv.getFeedcost()[1][0] - mkv.getBreedingsupplycost()[0] 
					- mkv.getHeiferpurchasecost()[0] - mkv.getVetcost()[0] - mkv.getLaborcost()[0] - mkv.getOthercost()[0] 
					- mkv.getFixedcost_labor()[0] - mkv.getFixedcost_other()[0]
					+ mkv.getCowsales()[0] - mkv.getLoss()[0]) / mkv.getCwt()[0];
				else mkv.getEaa_cwt()[y] = 0.0;
	
				mkv.getCash_cyr()[y] = mkv.getCash()[0] / years;	//not correct, because "repl" sales are included in t although cow is not present in t?
				if(mkv.getCwt()[0] > 0.0001) mkv.getCash_cwt()[y] = mkv.getCash()[0] / mkv.getCwt()[0];
				mkv.getMilkyield_cyr()[y] = mkv.getMilkyield()[0][0] / years;			//milk yield / cow / year
				if(mkv.getCwt()[0] > 0.0001) mkv.getMilkyield_cwt()[y] = mkv.getMilkyield()[0][0] / mkv.getCwt()[0];		//milk yield / cwt
				mkv.getMilksales_cyr()[y] = mkv.getMilksales()[0] / years;
				if(mkv.getCwt()[0] > 0.0001) mkv.getMilksales_cwt()[y] = mkv.getMilksales()[0] / mkv.getCwt()[0];
				mkv.getBreedingsupplycost_cyr()[y] = mkv.getBreedingsupplycost()[0] / years;
				if(mkv.getCwt()[0] > 0.0001) mkv.getBreedingsupplycost_cwt()[y] = mkv.getBreedingsupplycost()[0] / mkv.getCwt()[0];
				mkv.getFeedcost_cyr()[y] = (mkv.getFeedcost()[0][0] + mkv.getFeedcost()[1][0]) / years;
				if(mkv.getCwt()[0] > 0.0001) mkv.getFeedcost_cwt()[y] = (mkv.getFeedcost()[0][0] + mkv.getFeedcost()[1][0]) / mkv.getCwt()[0];
				mkv.getHeifercost_cyr()[y] = mkv.getHeiferpurchasecost()[0] / years;
				if(mkv.getCwt()[0] > 0.0001) mkv.getHeifercost_cwt()[y] = mkv.getHeiferpurchasecost()[0] / mkv.getCwt()[0];
				mkv.getCowsales_cyr()[y] = mkv.getCowsales()[0] / years;
				if(mkv.getCwt()[0] > 0.0001) mkv.getCowsales_cwt()[y] = mkv.getCowsales()[0] / mkv.getCwt()[0];
				mkv.getCalfsales_cyr()[y] = mkv.getCalfsales()[0] / years;
				if(mkv.getCwt()[0] > 0.0001) mkv.getCalfsales_cwt()[y] = mkv.getCalfsales()[0] / mkv.getCwt()[0];
				mkv.getLivestockcost_cyr()[y] = mkv.getHeifercost_cyr()[y] - mkv.getCowsales_cyr()[y];
				if(mkv.getCwt()[0] > 0.0001) mkv.getLivestockcost_cwt()[y] = mkv.getHeifercost_cwt()[y] - mkv.getCowsales_cwt()[y];
				mkv.getLaborcost_cyr()[y] = mkv.getLaborcost()[0] / years;
				if(mkv.getCwt()[0] > 0.0001) mkv.getLaborcost_cwt()[y] = mkv.getLaborcost()[0] / mkv.getCwt()[0];
				mkv.getLoss_cyr()[y] = mkv.getLoss()[0] / years;
				if(mkv.getCwt()[0] > 0.0001) mkv.getLoss_cwt()[y] = mkv.getLoss()[0] / mkv.getCwt()[0];
				mkv.getVetcost_cyr()[y] = mkv.getVetcost()[0] / years;
				if(mkv.getCwt()[0] > 0.0001) mkv.getVetcost_cwt()[y] = mkv.getVetcost()[0] / mkv.getCwt()[0];
				mkv.getOthercost_cyr()[y] = mkv.getOthercost()[0] / years;
				if(mkv.getCwt()[0] > 0.0001) mkv.getOthercost_cwt()[y] = mkv.getOthercost()[0] / mkv.getCwt()[0];
				mkv.getFixedcost_labor_cyr()[y] = mkv.getFixedcost_labor()[0] / years;
				if(mkv.getCwt()[0] > 0.0001) mkv.getFixedcost_labor_cwt()[y] = mkv.getFixedcost_labor()[0] / mkv.getCwt()[0];
				mkv.getFixedcost_other_cyr()[y] = mkv.getFixedcost_other()[0] / years;
				if(mkv.getCwt()[0] > 0.0001) mkv.getFixedcost_other_cwt()[y] = mkv.getFixedcost_other()[0] / mkv.getCwt()[0];
				mkv.getBreedingcost_cyr()[y] = (mkv.getBreedinglaborcost()[0] + mkv.getBreedingsupplycost()[0]) / years;
				if(mkv.getCwt()[0] > 0.0001) mkv.getBreedingcost_cwt()[y] = (mkv.getBreedinglaborcost()[0] + mkv.getBreedingsupplycost()[0]) / mkv.getCwt()[0];
	
				if(mkv.getPresent()[0][0] > 0.0) mkv.getCullrateall()[y] = (mkv.getRepl()[0][0] + mkv.getPind()[0][0]) / mkv.getPresent()[0][0] * Constants.SEA;
				if(mkv.getPresent()[0][0] > 0.0) mkv.getCullratevol()[y] = mkv.getRepl()[0][0] / mkv.getPresent()[0][0] * Constants.SEA;
				if(mkv.getPresent()[0][0] > 0.0) mkv.getCullrateinv()[y] = mkv.getPind()[0][0] / mkv.getPresent()[0][0] * Constants.SEA;
	
				if(mkv.getCullrateall()[y] > 0.0)	mkv.getHerdlife()[y] = 1/ mkv.getCullrateall()[y] * Constants.STAGE * Constants.SEA;
				if(mkv.getNewpreg()[0][0][0] > 0.0)	mkv.getPreglaborcost()[y] = mkv.getBreedinglaborcost()[0] / mkv.getNewpreg()[0][0][0];
				if(mkv.getNewpreg()[0][0][0] > 0.0)	mkv.getPregsupplycost()[y] = mkv.getBreedingsupplycost()[0] / mkv.getNewpreg()[0][0][0];
				if(mkv.getOpen()[0][0] > 0.0) 		mkv.getDaysopen()[0][y] = mkv.getDysopen()[0] / mkv.getOpen()[0][0] * Constants.STAGE;
	//&
				if(mkv.getNewpreg()[0][0][0] > 0.0)	mkv.getDaystoconc()[0][0][y][0] = mkv.getDaystoconc()[0][0][0][0] / mkv.getNewpreg()[0][0][0];
				if(h.getClvinterval_cnt()[0] > 0.0)	mkv.getClvinterval()[0][y][0] = mkv.getClvinterval()[0][0][0] / h.getClvinterval_cnt()[0];
				if(h.getDvwp_cnt()[0][0][0][0] > 0)	mkv.getDvwp()[0][0][y][0] = mkv.getDvwp()[0][0][0][0] / h.getDvwp_cnt()[0][0][0][0];
				if(h.getDfs_cnt()[0][0][0][0] > 0)	mkv.getDfs()[0][0][y][0] = mkv.getDfs()[0][0][0][0] / h.getDfs_cnt()[0][0][0][0];
				if(h.getDlbm_cnt()[0][0][0] > 0)	mkv.getDlbm()[0][0][y][0] = h.getDlbm()[0][0][0] / h.getDlbm_cnt()[0][0][0];
				
				if(mkv.getDysopencost_cnt()[0][0][0] > 0)	mkv.getDaysopen_cost()[y] = -mkv.getDysopencost()[0][0][0] / mkv.getDysopencost_cnt()[0][0][0];
				if(mkv.getDyspregcost_cnt()[0][0][0] > 0)	mkv.getDayspreg_cost()[y] = -mkv.getDyspregcost()[0][0][0] / mkv.getDyspregcost_cnt()[0][0][0];
				if(mkv.getPresent()[0][0] > 0.0) mkv.getOpendays()[y] = mkv.getOpen()[0][0] / mkv.getPresent()[0][0] * 365;
				if(mkv.getPresent()[0][0] - mkv.getMilking()[0] > 0.0) mkv.getFeedcostdry_cyr()[y] = mkv.getFeedcost()[0][0] / (mkv.getPresent()[0][0] - mkv.getMilking()[0]) * Constants.SEA;
				if(mkv.getMilking()[0] > 0.0) mkv.getFeedcostlac_cyr()[y] = mkv.getFeedcost()[1][0] / mkv.getMilking()[0] * Constants.SEA;
				if(mkv.getValue_cnt()[0] > 0.0) mkv.getCowvalue()[y] = mkv.getValue()[0] / mkv.getValue_cnt()[0];
				if(mkv.getMilking()[0] > 0.0) mkv.getDimy()[y] = mkv.getDim()[0] / mkv.getMilking()[0];
				if(mkv.getPind()[0][0]+mkv.getRepl()[0][0] > 0.0) mkv.getCullpricey()[y] = mkv.getCowsales()[0] / (mkv.getPind()[0][0] + mkv.getRepl()[0][0]);
				if(mkv.getAbortions_cnt()[0] > 0.0) mkv.getAbortrisk()[y] = mkv.getAbortions()[0] / mkv.getAbortions_cnt()[0] * 8;
				
	
				for(s = 0; s <= Constants.SEA; s++)
				{
					//pregnancy rate
					if(eligible_serv[s] > 0.0) mkv.getPregnancyrate0()[s][y] = eligible_preg[s] / eligible_serv[s];
					if(eligible_serv[s] + eligible_keep[s] > 0.0) mkv.getPregnancyrate1()[s][y] = eligible_preg[s] / (eligible_serv[s] + eligible_keep[s]);
	
					//fraction culled for "repl" and "pind" during the lactation, by calving season
	//				if(calved_clv[s] > 0.0) mkv.getfracrepl_clv[s][y] = repl_clv[s] / calved_clv[s];
	//				if(calved_clv[s] > 0.0) mkv.getfracpind_clv[s][y] = pind_clv[s] / calved_clv[s];
					//average dim when open cows got voluntary culled
	//				if(repl_open[s] > 0.0) mkv.getrepl_open_dim[s][y] = repl_open_dim[s] / repl_open[s];
				}
				//average pregnancy value (open, pregnant), average inseminationvalue
	
				if(pv0_cnt > 0.0) mkv.getPv_open()[y] = pv0 / pv0_cnt;
				if(pv10_cnt > 0.0) mkv.getPv_preg()[y] = pv10 / pv10_cnt;
	
				if(h.getInsval_cnt()[0][0] > 0.0) mkv.getInseminationvalue()[0][0][y] = h.getInsval()[0][0] / h.getInsval_cnt()[0][0];
	
				//statistics at culling open cows (annual averages)
				if(replopen_cnt > 0.0) mkv.getReplopen_cur()[0][0][y] = replopen_cur / replopen_cnt;
				if(replopen_cnt > 0.0) mkv.getReplopen_dim()[0][0][y] = replopen_dim / replopen_cnt;
				if(replopen_cnt > 0.0) mkv.getReplopen_mlk()[0][0][y] = replopen_mlk / replopen_cnt;
				if(replopen_cnt > 0.0) mkv.getReplopen_rev()[0][0][y] = replopen_rev / replopen_cnt;
				if(replopen_cnt > 0.0) mkv.getReplopen_iof()[0][0][y] = replopen_iof / replopen_cnt;
				if(replopen_cnt > 0.0) mkv.getReplopen_rpo()[0][0][y] = replopen_rpo / replopen_cnt;
	
	
	//*** reset variables so they can be used for the calculations for next year **********
				pvifa0 = 0.0;
				pvifa1 = 0.0;
				mkv.getAbortions()[0] = 0.0;
				mkv.getAbortions_cnt()[0] = 0.0;
				mkv.getBreedingsupplycost()[0] = 0.0;
				mkv.getBreedinglaborcost()[0] = 0.0;
				mkv.getCash()[0] = 0.0;
				mkv.getValue()[0] = 0.0;
				mkv.getValue_cnt()[0] = 0.0;
				mkv.getCullprice()[0] = 0.0;
				mkv.getCwt()[0] = 0.0;
				mkv.getMilkyield()[0][0] = 0.0;
				mkv.getMilksales()[0] = 0.0;
				mkv.getFeedcost()[0][0] = mkv.getFeedcost()[1][0] = 0.0;
				mkv.getHeiferpurchasecost()[0] = 0.0;
				mkv.getCowsales()[0] = 0.0;;
				mkv.getCalfsales()[0] = 0.0;
				mkv.getLaborcost()[0] = 0.0;
				mkv.getLoss()[0] = 0.0;
				mkv.getVetcost()[0] = 0.0;
				mkv.getOthercost()[0] = 0.0;
				mkv.getFixedcost_labor()[0] = 0.0;
				mkv.getFixedcost_other()[0] = 0.0;
				mkv.getClvinterval()[0][0][0] = 0.0;
				mkv.getCalved()[0] = 0.0;
	//&
				mkv.getClvinterval()[0][0][0] = 0.0;
				h.getClvinterval_cnt()[0] = 0.0;
				mkv.getDaystoconc()[0][0][0][0] = 0.0;
				mkv.getNewpreg()[0][0][0] = 0.0;
				mkv.getDvwp()[0][0][0][0] = 0.0;
				h.getDvwp_cnt()[0][0][0][0] = 0.0;
				mkv.getDfs()[0][0][0][0] = 0.0;
				h.getDfs_cnt()[0][0][0][0] = 0.0;
	
	
				mkv.getDysopencost()[0][0][0] = 0.0;
				mkv.getDysopencost_cnt()[0][0][0] = 0.0;
				mkv.getDyspregcost()[0][0][0] = 0.0;
				mkv.getDyspregcost_cnt()[0][0][0] = 0.0;
				mkv.getRepl()[0][0] = 0.0;
				mkv.getPind()[0][0] = 0.0;
				mkv.getPresent()[0][0] = 0.0;
				mkv.getMilking()[0] = 0.0;
				mkv.getOpen()[0][0] = 0.0;
				mkv.getDysopen()[0] = 0.0;
				mkv.getRpo()[0] = 0.0;
				mkv.getRpo_cnt()[0] = 0.0;
				mkv.getDim()[0] = 0.0;
	
	
				for(s = 0; s <= Constants.SEA; s++)
				{
					eligible_preg[s] = 0.0;
					eligible_keep[s] = 0.0;
					eligible_serv[s] = 0.0;
					repl_open_dim[s] = 0.0;
					repl_open[s] = 0.0;
				}
				repl_open[0] = 0.0;
				h.getInsval()[0][0] = 0.0;
				h.getInsval_cnt()[0][0] = 0.0;
				
				replopen_cnt = 0.0;
				replopen_dim = 0.0;
				replopen_cur = 0.0;
				replopen_mlk = 0.0;
				replopen_iof = 0.0;
				replopen_rev = 0.0;
				replopen_rpo = 0.0;
				pv0 = 0.0;
				pv0_cnt = 0.0;
				pv10 = 0.0;
				pv10_cnt = 0.0;
			}// end calculate statistics per year
			if((int)t % (int)((double)maxstages / 5) == 0 || t == 1 || t == maxstages)
				System.out.println((double)t/(double)maxstages*100);
			sea++; if(sea > Constants.SEA) sea = 1;
		}//end t
	
	//*** additional calculations **********
	
	
		sea = bio.getStart_sea();
		for(t = 1; t <= maxstages; t++)
		{
			for(mim = 1; mim <= Constants.MIM; mim++)
			{
				cs0 = calving_sea((short)(t + bio.getStart_sea() - mim - 1));		//start of mim: repl 
				cs1 = calving_sea((short)(t + bio.getStart_sea() - mim));			//end of mim: pind, calving
				ct0 = t - mim - 1;
				ct1 = t - mim;
				if(bio.getAnalysis() == 1) y = 1;	//t=1-12:y=1, t=13-24,y=2 		//steady state
				else
				{
					if(ct0 < 1) y = 0; else if(ct0 >=  1 && ct0 <= 12) y = 1; else if(ct0 >= 13 && ct0 <= 24) y = 2;
					else if(ct0 >= 25 && ct0 <= 36) y = 3; else if(ct0 >= 37) y = 4;
					if(ct1 < 1) y = 0; else if(ct1 >=  1 && ct1 <= 12) y = 1; else if(ct1 >= 13 && ct1 <= 24) y = 2;
					else if(ct1 >= 25 && ct1 <= 36) y = 3; else if(ct1 >= 37) y = 4;
				}
	//&		//System.out.println("\nMARKOV_CHAIN_STATS t:%d sea:%2d mim:%2d  cs0:%2d cs1:%2d ct0:%2d ct1:%2d y:%d",t,sea,mim,cs0,cs1,ct0,ct1,y); 
				for(lact = 1; lact <= Constants.LACT; lact++)
				{
					mkv.getCalvinginterval_clv()[lact][0] += mkv.getClvinterval()[lact][mim][t] * (double)mim;
					mkv.getCalvinginterval_clv()[0][cs1] += mkv.getClvinterval()[lact][mim][t] * (double)mim;
					calvinginterval_clv_cnt[lact][0] += mkv.getClvinterval()[lact][mim][t];
					calvinginterval_clv_cnt[0][cs1] += mkv.getClvinterval()[lact][mim][t];
	
					for(curr = 1; curr <= Constants.CURR; curr++) 
					{
						mkv.getDaystoconception_clv()[lact][0] += mkv.getDaystoconc()[curr][lact][mim][t] * (double)mim * Constants.STAGE;
						mkv.getDaystoconception_clv()[0][cs1] += mkv.getDaystoconc()[curr][lact][mim][t] * (double)mim * Constants.STAGE;
						daystoconception_clv_cnt[lact][0] += mkv.getDaystoconc()[curr][lact][mim][t];
						daystoconception_clv_cnt[0][cs1] += mkv.getDaystoconc()[curr][lact][mim][t];
	
						mkv.getDaysvwp_clv()[lact][0] += mkv.getDvwp()[curr][lact][mim][t] * (double)mim * Constants.STAGE;
						mkv.getDaysvwp_clv()[0][cs1] += mkv.getDvwp()[curr][lact][mim][t] * (double)mim * Constants.STAGE;
						daysvwp_clv_cnt[lact][0] += mkv.getDvwp()[curr][lact][mim][t];
						daysvwp_clv_cnt[0][cs1] += mkv.getDvwp()[curr][lact][mim][t];
	
						mkv.getDaysfs_clv()[lact][0] += mkv.getDfs()[curr][lact][mim][t] * (double)mim * Constants.STAGE;
						mkv.getDaysfs_clv()[0][cs1] += mkv.getDfs()[curr][lact][mim][t] * (double)mim * Constants.STAGE;
						daysfs_clv_cnt[lact][0] += mkv.getDfs()[curr][lact][mim][t];
						daysfs_clv_cnt[0][cs1] += mkv.getDfs()[curr][lact][mim][t];
	
						//mkv.getdayslbm_clv[0][cs1] not calculated because quite complicated
					}
				}
				pind_clv[cs1][y] += mkv.getPind_mim()[mim][t] * (double)mim;					//end of stage
				pind_clv_cnt[cs1][y] += mkv.getPind_mim()[mim][t];
				repl_clv[cs0][y] += mkv.getRepl_mim()[mim][t] * (double)(mim - 1);				//begin of stage
				repl_clv_cnt[cs0][y] += mkv.getRepl_mim()[mim][t];
			}
	
	//stop(0);
	
			sea++; if(sea > Constants.SEA) sea = 1;
		}//end t
	
	
		for(mim = 1; mim <= Constants.MIM; mim++)
		{
			for(t = 1; t <= maxstages; t++)
			{
				mkv.getEligible_bred()[mim][0] += mkv.getEligible_bred()[mim][t];
				mkv.getEligible_keep()[mim][0] += mkv.getEligible_keep()[mim][t];
				mkv.getEligible_serv()[mim][0] += mkv.getEligible_serv()[mim][t];
				mkv.getEligible_preg()[mim][0] += mkv.getEligible_preg()[mim][t];
			}
			mkv.getEligible_bred()[mim][0] = mkv.getEligible_bred()[mim][0] / (double)maxstages;
			mkv.getEligible_keep()[mim][0] = mkv.getEligible_keep()[mim][0] / (double)maxstages;
			mkv.getEligible_serv()[mim][0] = mkv.getEligible_serv()[mim][0] / (double)maxstages;
			mkv.getEligible_preg()[mim][0] = mkv.getEligible_preg()[mim][0] / (double)maxstages;
		}
	
	
		//System.out.println("lact : "+lact);
		//System.out.println("h.getDlbm_cnt()[0][lact][0] > 0.0:  "+(h.getDlbm_cnt()[0][lact][0] > 0.0));
		//System.out.println("mkv.getDayslbm_clv()[lact][0]:  "+mkv.getDayslbm_clv()[lact][0]);
		//System.out.println("h.getDlbm()[0][lact][0]:  "+h.getDlbm()[0][lact][0]);
		//System.out.println("h.getDlbm_cnt()[0][lact][0]:  "+h.getDlbm_cnt()[0][lact][0]);
		
		/**
		 * 
		 * I suppose this should not be present at all ....
		 * 
		 * if(h.getDlbm_cnt()[0][lact][0] > 0.0) mkv.getDayslbm_clv()[lact][0] = h.getDlbm()[0][lact][0] / h.getDlbm_cnt()[0][lact][0];
		 * 
		 */
		//	System.out.println("\n lact:%d  mkv.getDlbm():%f",lact,mkv.getDlbm()[0][lact][0][0]);
	
	
	
	
	
	
	
	
	//&
		for(lact = 1; lact <= Constants.LACT; lact++) if(daystoconception_clv_cnt[lact][0] > 0.0) 
			mkv.getDaystoconception_clv()[lact][0] = mkv.getDaystoconception_clv()[lact][0] / daystoconception_clv_cnt[lact][0];
		for(s = 1; s <= Constants.SEA; s++) if(daystoconception_clv_cnt[0][s] > 0.0)
			mkv.getDaystoconception_clv()[0][s] = mkv.getDaystoconception_clv()[0][s] / daystoconception_clv_cnt[0][s];
	
		for(lact = 1; lact <= Constants.LACT; lact++) if(calvinginterval_clv_cnt[lact][0] > 0.0) 
			mkv.getCalvinginterval_clv()[lact][0] = mkv.getCalvinginterval_clv()[lact][0] / calvinginterval_clv_cnt[lact][0];
		for(s = 1; s <= Constants.SEA; s++) if(calvinginterval_clv_cnt[0][s] > 0.0)
			mkv.getCalvinginterval_clv()[0][s] = mkv.getCalvinginterval_clv()[0][s] / calvinginterval_clv_cnt[0][s];
	
		for(lact = 1; lact <= Constants.LACT; lact++) if(daysvwp_clv_cnt[lact][0] > 0.0) 
			mkv.getDaysvwp_clv()[lact][0] = mkv.getDaysvwp_clv()[lact][0] / daysvwp_clv_cnt[lact][0];
		for(s = 1; s <= Constants.SEA; s++) if(daysvwp_clv_cnt[0][s] > 0.0)
			mkv.getDaysvwp_clv()[0][s] = mkv.getDaysvwp_clv()[0][s] / daysvwp_clv_cnt[0][s];
	
		for(lact = 1; lact <= Constants.LACT; lact++) if(daysfs_clv_cnt[lact][0] > 0.0) 
			mkv.getDaysfs_clv()[lact][0] = mkv.getDaysfs_clv()[lact][0] / daysfs_clv_cnt[lact][0];
		for(s = 1; s <= Constants.SEA; s++) if(daysfs_clv_cnt[0][s] > 0.0)
			mkv.getDaysfs_clv()[0][s] = mkv.getDaysfs_clv()[0][s] / daysfs_clv_cnt[0][s];
	
		for(lact = 1; lact <= Constants.LACT; lact++)
		{
			for(curr = 1; curr <= Constants.CURR; curr++)
			{
				for(t = 1; t <= maxstages; t++)
				{
					h.getDlbm()[0][lact][0] += h.getDlbm()[curr][lact][t];
					h.getDlbm_cnt()[0][lact][0] += h.getDlbm_cnt()[curr][lact][t];
				}
			}
			if(h.getDlbm_cnt()[0][lact][0] > 0.0) mkv.getDayslbm_clv()[lact][0] = h.getDlbm()[0][lact][0] / h.getDlbm_cnt()[0][lact][0];
		//	System.out.println("\n lact:%d  mkv.getDlbm():%f",lact,mkv.getDlbm()[0][lact][0][0]);
		}
		for(curr = 1; curr <= Constants.CURR; curr++)
		{
			for(lact = 1; lact <= Constants.LACT; lact++)
			{
				for(t = 1; t <= maxstages; t++)
				{
					h.getDlbm()[curr][0][0] += h.getDlbm()[curr][lact][t];
					h.getDlbm_cnt()[curr][0][0] += h.getDlbm_cnt()[curr][lact][t];
				//	System.out.println("\n curr:%d lact:%d t:%d  h.getDlbm():%f  h.getDlbm_cnt():%f",curr,lact,t,h.getDlbm()[curr][lact][t],h.getDlbm_cnt()[curr][lact][t]);
				//	System.out.println(" %f",h.getDlbm()[curr][0][0]/h.getDlbm_cnt()[curr][0][0]);
				}
			//	stop(0);
			}
			if(h.getDlbm_cnt()[curr][0][0] > 0.0) mkv.getDlbm()[curr][0][0][0] = h.getDlbm()[curr][0][0] / h.getDlbm_cnt()[curr][0][0];
		//	System.out.println("\n curr:%d  mkv.getDlbm():%f",curr,mkv.getDlbm()[curr][0][0][0]);
		}
	//	stop(0);
	
	
	
	
		for(s = 1; s <= Constants.SEA; s++)
			for(y = 1; y <= Constants.YEAR; y++)
			{
				if(pind_clv_cnt[s][y] > 0.0) mkv.getPind_clv()[s][y] += pind_clv[s][y] / pind_clv_cnt[s][y];
				if(repl_clv_cnt[s][y] > 0.0) mkv.getRepl_clv()[s][y] += repl_clv[s][y] / repl_clv_cnt[s][y];
			}
	
	//	//fraction open and pregnant cows by mim in t = maxstages + 1 (=typically start season)
	//replaced by mkv.getStatus()[mim][] in markov_chain.
	//	sum = 0.0;
	//	for(curr = 1; curr <= Constants.CURR; curr++)
	//		for(lact = 1; lact <= bio.getMaxlact(); lact++)
	//			for(mim = 1; mim <= Constants.MIM; mim++)
	//				for(mpreg = 0; mpreg <= Constants.MPREG; mpreg++)
	//				{
	//					sum += mkv.getcow[curr][lact][mim][mpreg][2];
	//					if(mpreg == 0) mkv.getStatus()[mim][0] += mkv.getcow[curr][lact][mim][mpreg][2];
	//					else mkv.getStatus()[mim][1] += mkv.getcow[curr][lact][mim][mpreg][2];
	//					
	//				//	System.out.println("\n curr:%d lact:%d mim:%d mpreg:%d cow:%f",curr,lact,mim,mpreg,mkv.getcow[curr][lact][mim][mpreg][2]);
	//				//	System.out.println(" open:%f preg:%f",mkv.getStatus()[mim][0],mkv.getStatus()[mim][1]);
	//				//	if(mim <= 2 && mpreg == 2 && mkv.getStatus()[2][1] > 0.0) stop(0);
	//				}
	//	mkv.getStatus()[1][0] += mkv.getheiferpool[1];
	//mkv.getcow[] does not include heifer. Adding heiferpool[1] assumes that the herd is in 
	//steady state and season at t=1 == season at maxstage + 1
	//This assumption only affects a graph in newdairy1.xls and no herd statistics.
	
	
		for(mim = 1; mim <= Constants.MIM; mim++)
		{
			mkv.getStatus()[mim][0] = mkv.getStatus()[mim][0] / (double)maxstages;
			mkv.getStatus()[mim][1] = mkv.getStatus()[mim][1] / (double)maxstages;
		}
	
	
		//risk of culling per month in lactation
		for(mim = 1; mim <= Constants.MIM; mim++)
		{
			h.getAtrisk()[mim] = mkv.getHz_preg()[mim][0] + mkv.getHz_open()[mim][0];
			h.getCulled()[mim] = mkv.getHz_open()[mim][1] + mkv.getHz_open()[mim][2] + mkv.getHz_preg()[mim][1] + mkv.getHz_preg()[mim][2];
			h.getCulled()[0]  += mkv.getHz_open()[mim][1] + mkv.getHz_open()[mim][2] + mkv.getHz_preg()[mim][1] + mkv.getHz_preg()[mim][2];
			if(h.getAtrisk()[mim] > 0.0) mkv.getHz_cull()[mim] = h.getCulled()[mim] / h.getAtrisk()[mim];
	
		//	System.out.println("\nopen %f > %f ",mkv.getHz_open()[mim][0],mkv.getHz_open()[mim][1]+mkv.getHz_open()[mim][2]),
		//	System.out.println("  preg %f > %f ",mkv.getHz_preg()[mim][0],mkv.getHz_preg()[mim][1]+mkv.getHz_preg()[mim][2]),
		//	System.out.println("   all %f > %f  risk:%f",h.getAtrisk()[mim],h.getCulled()[mim],mkv.getHz_cull()[mim]);
	
			if(mkv.getHz_open()[mim][0] > 0.0) mkv.getHz_open()[mim][1] = mkv.getHz_open()[mim][1] / mkv.getHz_open()[mim][0];	//repl open cows
			if(mkv.getHz_open()[mim][0] > 0.0) mkv.getHz_open()[mim][2] = mkv.getHz_open()[mim][2] / mkv.getHz_open()[mim][0];	//pind open cows
			if(mkv.getHz_preg()[mim][0] > 0.0) mkv.getHz_preg()[mim][1] = mkv.getHz_preg()[mim][1] / mkv.getHz_preg()[mim][0];	//repl pregnant cows
			if(mkv.getHz_preg()[mim][0] > 0.0) mkv.getHz_preg()[mim][2] = mkv.getHz_preg()[mim][2] / mkv.getHz_preg()[mim][0];	//pind pregnant cows
	
	//if(mkv.getHz_preg()[mim][1] > 0.0)
	//{
	//	System.out.println("\nhz_preg[%d][1]=%.10f > 0",mim,mkv.getHz_preg()[mim][1]);
	//	stop(0);
	//}
	
		}
		for(mim = 1; mim <= Constants.MIM; mim++) 
		{
			if(h.getCulled()[0] > 0.01)
				mkv.getCulled()[mim] = h.getCulled()[mim] / h.getCulled()[0];
			//System.out.println("\n culled:%d %f %f %f ",mim,mkv.getCulled()[mim],h.getCulled()[mim],h.getCulled()[0]);
		}
	//&
		//days to conception per curr
		//days to vwp per curr
		for(curr = 1; curr <= Constants.CURR; curr++)
		{
			if(mkv.getNewpreg()[curr][0][0]   > 0.0) mkv.getDaystoconc()[curr][0][0][0] = mkv.getDaystoconc()[curr][0][0][0] / mkv.getNewpreg()[curr][0][0];
			if(h.getDvwp_cnt()[curr][0][0][0] > 0.0) mkv.getDvwp()[curr][0][0][0] = mkv.getDvwp()[curr][0][0][0] / h.getDvwp_cnt()[curr][0][0][0];
			if(h.getDfs_cnt()[curr][0][0][0] > 0.0) mkv.getDfs()[curr][0][0][0] = mkv.getDfs()[curr][0][0][0] / h.getDfs_cnt()[curr][0][0][0];
	
		//	System.out.println("\n curr:%d %f %f",curr,mkv.getDaystoconc()[curr][0][0] ,mkv.getdaystoconception[curr][0]);
		}
		//stop(0);
	
		for(curr = 1; curr <= Constants.CURR; curr++)
			if(mkv.getDysopencost_cnt()[curr][0][0] > 0.0) mkv.getDysopencost()[curr][0][0] = -mkv.getDysopencost()[curr][0][0] / mkv.getDysopencost_cnt()[curr][0][0];
		for(mim = 1; mim <= Constants.MIM; mim++)
			if(mkv.getDysopencost_cnt()[0][mim][0] > 0.0) mkv.getDysopencost()[0][mim][0] = -mkv.getDysopencost()[0][mim][0] / mkv.getDysopencost_cnt()[0][mim][0];
	
		for(curr = 1; curr <= Constants.CURR; curr++)
			if(mkv.getDyspregcost_cnt()[curr][0][0] > 0.0) mkv.getDyspregcost()[curr][0][0] = -mkv.getDyspregcost()[curr][0][0] / mkv.getDyspregcost_cnt()[curr][0][0];
		for(mim = 1; mim <= Constants.MIM; mim++)
			if(mkv.getDyspregcost_cnt()[0][mim][0] > 0.0) mkv.getDyspregcost()[0][mim][0] = -mkv.getDyspregcost()[0][mim][0] / mkv.getDyspregcost_cnt()[0][mim][0];
	
		//statistics at voluntary culling
		if(mkv.getReplopen_cnt()[0][0][0] > 0.0) mkv.getReplopen_dim()[0][0][0] = mkv.getReplopen_dim()[0][0][0] / mkv.getReplopen_cnt()[0][0][0];
		if(mkv.getReplopen_cnt()[0][0][0] > 0.0) mkv.getReplopen_cur()[0][0][0] = mkv.getReplopen_cur()[0][0][0] / mkv.getReplopen_cnt()[0][0][0];
		if(mkv.getReplopen_cnt()[0][0][0] > 0.0) mkv.getReplopen_rev()[0][0][0] = mkv.getReplopen_rev()[0][0][0] / mkv.getReplopen_cnt()[0][0][0];
		if(mkv.getReplopen_cnt()[0][0][0] > 0.0) mkv.getReplopen_iof()[0][0][0] = mkv.getReplopen_iof()[0][0][0] / mkv.getReplopen_cnt()[0][0][0];
		if(mkv.getReplopen_cnt()[0][0][0] > 0.0) mkv.getReplopen_mlk()[0][0][0] = mkv.getReplopen_mlk()[0][0][0] / mkv.getReplopen_cnt()[0][0][0];
		if(mkv.getReplopen_cnt()[0][0][0] > 0.0) mkv.getReplopen_rpo()[0][0][0] = mkv.getReplopen_rpo()[0][0][0] / mkv.getReplopen_cnt()[0][0][0];
		for(mim = 1; mim <= Constants.MIM; mim++)
		{
			if(mkv.getReplopen_cnt()[mim][0][0] > 0.0) mkv.getReplopen_dim()[mim][0][0] = mkv.getReplopen_dim()[mim][0][0] / mkv.getReplopen_cnt()[mim][0][0];
			if(mkv.getReplopen_cnt()[mim][0][0] > 0.0) mkv.getReplopen_cur()[mim][0][0] = mkv.getReplopen_cur()[mim][0][0] / mkv.getReplopen_cnt()[mim][0][0];
			if(mkv.getReplopen_cnt()[mim][0][0] > 0.0) mkv.getReplopen_rev()[mim][0][0] = mkv.getReplopen_rev()[mim][0][0] / mkv.getReplopen_cnt()[mim][0][0];
			if(mkv.getReplopen_cnt()[mim][0][0] > 0.0) mkv.getReplopen_iof()[mim][0][0] = mkv.getReplopen_iof()[mim][0][0] / mkv.getReplopen_cnt()[mim][0][0];
			if(mkv.getReplopen_cnt()[mim][0][0] > 0.0) mkv.getReplopen_mlk()[mim][0][0] = mkv.getReplopen_mlk()[mim][0][0] / mkv.getReplopen_cnt()[mim][0][0];
			if(mkv.getReplopen_cnt()[mim][0][0] > 0.0) mkv.getReplopen_rpo()[mim][0][0] = mkv.getReplopen_rpo()[mim][0][0] / mkv.getReplopen_cnt()[mim][0][0];
		}
		for(t = 1; t <= Constants.MAXM; t++)
		{
			if(mkv.getReplopen_cnt()[0][t][0] > 0.0) mkv.getReplopen_dim()[0][t][0] = mkv.getReplopen_dim()[0][t][0] / mkv.getReplopen_cnt()[0][t][0];
			if(mkv.getReplopen_cnt()[0][t][0] > 0.0) mkv.getReplopen_cur()[0][t][0] = mkv.getReplopen_cur()[0][t][0] / mkv.getReplopen_cnt()[0][t][0];
			if(mkv.getReplopen_cnt()[0][t][0] > 0.0) mkv.getReplopen_rev()[0][t][0] = mkv.getReplopen_rev()[0][t][0] / mkv.getReplopen_cnt()[0][t][0];
			if(mkv.getReplopen_cnt()[0][t][0] > 0.0) mkv.getReplopen_iof()[0][t][0] = mkv.getReplopen_iof()[0][t][0] / mkv.getReplopen_cnt()[0][t][0];
			if(mkv.getReplopen_cnt()[0][t][0] > 0.0) mkv.getReplopen_mlk()[0][t][0] = mkv.getReplopen_mlk()[0][t][0] / mkv.getReplopen_cnt()[0][t][0];
			if(mkv.getReplopen_cnt()[0][t][0] > 0.0) mkv.getReplopen_rpo()[0][t][0] = mkv.getReplopen_rpo()[0][t][0] / mkv.getReplopen_cnt()[0][t][0];
		}
	
	
		//inseminationvalue by mim
		for(mim = 1; mim <= Constants.MIM; mim++)
		{
			if(h.getInsval_cnt()[mim][0] > 0.0)     mkv.getInseminationvalue()[mim][0][0] = h.getInsval()[mim][0] / h.getInsval_cnt()[mim][0];
		}
		//pinseminationvalue by t
		for(t = 1; t <= maxstages; t++)
		{
			if(h.getInsval_cnt()[0][t] > 0.0)     mkv.getInseminationvalue()[0][t][0] = h.getInsval()[0][t] / h.getInsval_cnt()[0][t];
		}
	
	
		//pregnancyvalue by curr (open cows: mpreg = 0,  pregnant cows: mpreg = 10)
		for(curr = 1; curr <= Constants.CURR; curr++) for(mpreg = 0; mpreg <= Constants.MPREG + 1; mpreg++) if(mkv.getPv_curr_cnt()[curr][mpreg] > 0)
			mkv.getPv_curr()[curr][mpreg] = mkv.getPv_curr()[curr][mpreg] / mkv.getPv_curr_cnt()[curr][mpreg];
	
		//pregnancyvalue by lact (open cows: mpreg = 0,  pregnant cows: mpreg = 10)
		for(lact = 1; lact <= bio.getMaxlact(); lact++) for(mpreg = 0; mpreg <= Constants.MPREG + 1; mpreg++) if(mkv.getPv_lact_cnt()[lact][mpreg] > 0)
			mkv.getPv_lact()[lact][mpreg] = mkv.getPv_lact()[lact][mpreg] / mkv.getPv_lact_cnt()[lact][mpreg];
	
		//pregnancyvalue by mim (open cows: mpreg = 0,  pregnant cows: mpreg = 10)
		for(mim = 1; mim <= Constants.MIM; mim++) for(mpreg = 0; mpreg <= Constants.MPREG + 1; mpreg++) if(mkv.getPv_mim_cnt()[mim][mpreg] > 0)
			mkv.getPv_mim()[mim][mpreg] = mkv.getPv_mim()[mim][mpreg] / mkv.getPv_mim_cnt()[mim][mpreg];
	
		//pregnancyvalue by mpreg (open cows: mpreg = 0,  pregnant cows: mpreg = 10)
		for(mpreg = 0; mpreg <= Constants.MPREG + 1; mpreg++) if(mkv.getPv_mpreg_cnt()[mpreg] > 0)
			mkv.getPv_mpreg()[mpreg] = mkv.getPv_mpreg()[mpreg] / mkv.getPv_mpreg_cnt()[mpreg];
	
		//pregnancyvalue by t (open cows: mpreg = 0,  pregnant cows: mpreg = 10)
		for(t = 1; t <= maxstages; t++) for(mpreg = 0; mpreg <= Constants.MPREG + 1; mpreg++) if(mkv.getPv_t_cnt()[t][mpreg] > 0)
			mkv.getPv_t()[t][mpreg] = mkv.getPv_t()[t][mpreg] / mkv.getPv_t_cnt()[t][mpreg];
	
		//pregnancyvalue by curr and lact (open cows: mpreg = 0,  pregnant cows: mpreg = 10)
		for(curr = 1; curr <= Constants.CURR; curr++) for(lact = 1; lact <= bio.getMaxlact(); lact++) for(mpreg = 0; mpreg <= Constants.MPREG + 1; mpreg++) if(mkv.getPv_curr_lact_cnt()[curr][lact][mpreg] > 0)
			mkv.getPv_curr_lact()[curr][lact][mpreg] = mkv.getPv_curr_lact()[curr][lact][mpreg] / mkv.getPv_curr_lact_cnt()[curr][lact][mpreg];
	
	
		//dry-off yield
		for(t = 1; t <= maxstages; t++)
			if(mkv.getDryoffyield_cnt()[t] > 0.0) 
				mkv.getDryoffyield()[t] = mkv.getDryoffyield()[t] / mkv.getDryoffyield_cnt()[t];
											
		//lactation curve + distribution of cows by mim
		for(mim = 1; mim <= Constants.MIM; mim++)
		{
			for(t = 1; t <= maxstages; t++)
			{
				mkv.getMilkyield()[mim][0] += mkv.getMilkyield()[mim][t];
				mkv.getMilkyield_cnt()[mim][0] += mkv.getMilkyield_cnt()[mim][t];
				mkv.getMilkyield_cnt()[0][0] += mkv.getMilkyield_cnt()[mim][t];
			//	System.out.println("\n mim:%d  t:%d  %f %f %f",mim,t,mkv.getMilkyield()[mim][t],mkv.getMilkyield_cnt()[mim][t],mkv.getMilkyield_cnt()[0][0]);
			//	stop(0);
			}
			if(mkv.getMilkyield_cnt()[mim][0] > 0)
				mkv.getMilkyield()[mim][0] = mkv.getMilkyield()[mim][0] / mkv.getMilkyield_cnt()[mim][0] / Constants.STAGE;
		}
		for(mim = 1; mim <= Constants.MIM; mim++)
			if(mkv.getMilkyield_cnt()[0][0] > 0)
				mkv.getMilkyield_cnt()[mim][0] = mkv.getMilkyield_cnt()[mim][0] / mkv.getMilkyield_cnt()[0][0];
	
	
		//%cows per lactation number
		for(lact = 1; lact <= bio.getMaxlact(); lact++)
		{
			for(t = 1; t <= maxstages; t++)
			{
				mkv.getLact()[lact][0] += mkv.getLact()[lact][t];
				mkv.getLact()[0][0] += mkv.getLact()[lact][t];
			}
		}
		for(lact = 1; lact <= bio.getMaxlact(); lact++)
			if(mkv.getLact()[0][0] > 0) mkv.getLact()[lact][0] = mkv.getLact()[lact][0] / mkv.getLact()[0][0];
	
	
		//pind_rpo (ISVEE Cairns)
		for(curr = 1; curr <= Constants.CURR; curr++) for(lact = 1; lact <= bio.getMaxlact(); lact++) for(mim = 1; mim <= Constants.MIM; mim++) for(mpreg = 0; mpreg <= Constants.MPREG; mpreg++)	
			if(mkv.getPind_rpo_cnt()[curr][lact][mim][mpreg] > 0.0)
			{	
				//curr
				mkv.getPind_rpo()[curr][0][0][0] += mkv.getPind_rpo()[curr][lact][mim][mpreg];
				mkv.getPind_rpo_cnt()[curr][0][0][0] += mkv.getPind_rpo_cnt()[curr][lact][mim][mpreg];
	
				//lact
				mkv.getPind_rpo()[0][lact][0][0] += mkv.getPind_rpo()[curr][lact][mim][mpreg];
				mkv.getPind_rpo_cnt()[0][lact][0][0] += mkv.getPind_rpo_cnt()[curr][lact][mim][mpreg];
					
				//mim (open cows)
				if(mpreg == 0) mkv.getPind_rpo()[0][0][mim][0] += mkv.getPind_rpo()[curr][lact][mim][0];
				if(mpreg == 0) mkv.getPind_rpo_cnt()[0][0][mim][0] += mkv.getPind_rpo_cnt()[curr][lact][mim][0];
			
				//mpreg
				mkv.getPind_rpo()[0][0][0][mpreg] += mkv.getPind_rpo()[curr][lact][mim][mpreg];
				mkv.getPind_rpo_cnt()[0][0][0][mpreg] += mkv.getPind_rpo_cnt()[curr][lact][mim][mpreg];
	
			//	System.out.println("\n %d %d %d %d %f %f %f %f",curr,lact,mim,mpreg,
			//		mkv.getPind_rpo()[curr][0][0][0],mkv.getPind_rpo()[0][lact][0][0],mkv.getPind_rpo()[0][0][mim][0],mkv.getPind_rpo()[0][0][0][mpreg]);
			//	if(mkv.getPind_rpo()[curr][lact][mim][mpreg] > 0.0) stop(1);
			}
		for(curr = 1; curr <= Constants.CURR; curr++) if(mkv.getPind_rpo_cnt()[curr][0][0][0] > 0.0)
			mkv.getPind_rpo()[curr][0][0][0] = mkv.getPind_rpo()[curr][0][0][0] / mkv.getPind_rpo_cnt()[curr][0][0][0];
		
		for(lact = 1; lact <= bio.getMaxlact(); lact++) if(mkv.getPind_rpo_cnt()[0][lact][0][0] > 0.0)
			mkv.getPind_rpo()[0][lact][0][0] = mkv.getPind_rpo()[0][lact][0][0] / mkv.getPind_rpo_cnt()[0][lact][0][0];
		
		for(mim = 1; mim <= Constants.MIM; mim++) if(mkv.getPind_rpo_cnt()[0][0][mim][0] > 0.0)
			mkv.getPind_rpo()[0][0][mim][0] = mkv.getPind_rpo()[0][0][mim][0] / mkv.getPind_rpo_cnt()[0][0][mim][0];
	
		for(mpreg = 0; mpreg <= Constants.MPREG; mpreg++) if(mkv.getPind_rpo_cnt()[0][0][0][mpreg] > 0.0)
			mkv.getPind_rpo()[0][0][0][mpreg] = mkv.getPind_rpo()[0][0][0][mpreg] / mkv.getPind_rpo_cnt()[0][0][0][mpreg];
	
	
		for(mim = 1; mim <= Constants.MIM; mim++)
		{
			sum = 0.0;
			for(t = 1; t <= maxstages; t++)
			{
				mkv.getDmi()[mim][0] += mkv.getDmi()[mim][t];
				sum += mkv.getDmi_cnt()[mim][t];
			}
			if(sum > 0.0) mkv.getDmi()[mim][0] = mkv.getDmi()[mim][0] / sum;
		}
	
	
	//	avg = 0.0; System.out.println("\n discounted\n");
	//	for(t = 1; t <= Constants.SEA; t++)
	//	{
	//		System.out.println("%.2f + ",mkv.getCash()[t]);	
	//		avg += mkv.getCash()[t];
	//	}
	//	System.out.println(" = cash:%f",avg);	//should equal fut_cash in DP_POLICY()
	//	avg = 0.0; System.out.println("\n not discounted\n");
	//	for(t = 1; t <= Constants.SEA; t++)
	//	{
	//		System.out.println("%.2f + ",mkv.getCash()[t] * pow(1+bio.getInterest()/Constants.SEA,(double)t-1));	
	//		avg += mkv.getCash()[t] * pow(1+bio.getInterest()/Constants.SEA,(double)t-1);
	//	}
	//	System.out.println(" = cash:%f",avg);	//should equal fut_cash in DP_POLICY()
	//	stop(0);
	
		System.out.println("\n MARKOV_CHAIN_STATS done ...");
		Utility.free(h);
		h = null;
		return;
	}/*** end MARKOV_CHAIN_STATS ***/
	
	
	/******************************************************************************
	PREGNANCY_VALUE: Determines the pregnancy value of a cow.
	- see Dairycomp305 for CowValue and PregValue definitions (search Google)
	- http://www.puyallup.wsu.edu/dairy/shortcourse/readArticle.asp?intArticleID=74
	
	Pregnancy Value = FV of pregnant cow compared to FV open cow = prg
	- If cow is open: prg = rpo(pregnant) - rpo(open) = extra value of cow if she conceived now:
	  rpo.getprg[curr][lact][mim][0][sea] = rpo.getcow[curr][lact][mim][1][sea] - rpo.getcow[curr][lact][mim][0][sea];
	
	- If cow is pregnant: prg = rpo(pregnant) - rpo(open) = loss in value if cow aborted now
	  MPREG = 1-5: same lactation continues:
	  rpo.getprg[curr][lact][mim][mpreg][sea] = rpo.getcow[curr][lact][mim][mpreg][sea] - rpo.getcow[curr][lact][mim][0][sea];
	
	  MPREG = 6-9: next lactation:
	  rpo.getprg[curr][lact][mim][mpreg][sea] = rpo.getcow[curr][lact][mim][mpreg][sea] - rpo.getcow[curr][lact+1][1][0][sea];
	
	Note: ins.getcow[] is Insemination Value = value of breeding the open cow vs not breeding her (at begin of month)
	      see Dairycomp305 for CowValue and PregValue definitions (search Google).
	
	(question is, would this cow continue her current lactation or start the next lactation? may depend on mpreg: see below)
	Note: every cow that starts a new lactation has a calf to sell (DP_POLICY) -.get add "abort" state?
	Note: how is milk yield affected of cows that start a new lactation because of culling?
	
	13. Cows Aborting or Calving Prematurely (http://www.dhia.org/udcp.pdf)
	Cows freshening 30 or more days prior to the expected calving date, whether in milk or dry, will be coded as abnormal
	(abortion). When a breeding date is available, a cow calving less than 30 days prior to the expected calving date will
	be considered a normal calving. If a cow aborts while in milk and has carried a calf less than 152 days (5 months), 
	her current record will continue without interruption. If a breeding date is not available, and the cow aborts while
	in milk for less than 200 days (6.5 months), her current record will continue without interruption. Except for the
	specific situations above, the current record will end and a new lactation will begin.
	******************************************************************************/
	public void pregnancy_value(BioRecord bio, RpoRecord rpo)
	{
		int curr = 0, lact = 0, mim = 0, mpreg = 0, sea = 0, lact3 = 0;
	
	//Note: build abortion in model (mpreg < MPREG returns to mpreg = 0)?
	//Calving with MPREG == 8 may give viable calf -.get add to prg.getvalue?
	//mpreg > 5 -.get next lactation, same curr level?  (PregVal (DC305) uses same lactation)
	//mpreg > 5, lact == LACT -.get next lactation not available
	
		System.out.println("\n PREGNANCY_VALUE started ...");
		for(curr = 1; curr <= Constants.CURR; curr++)
		{
			for(lact = 1; lact <= bio.getMaxlact(); lact++)
			{
				if(lact <= 3) lact3 = lact; else lact3 = 3;
				for(mim = 1; mim <= Constants.MIM; mim++)
				{
					for(sea = 1; sea <= Constants.SEA; sea++)
					{
						for(mpreg = 0; mpreg <= Constants.MPREG; mpreg++) if(mpreg == 0 || (mim - mpreg >= bio.getVwp()[lact3] && mim - mpreg <= bio.getLbm()))
						{
						//	System.out.println("\n curr:%d lact:%d mim:%d sea:%d mpreg:%d rpo:%.0f",curr,lact,mim,sea,mpreg,rpo.getcow[curr][lact][mim][mpreg][sea]);
						//	if(mpreg <= 5)			//open cows + early abortion, same lactation
							if(mpreg <= Constants.MPREG)
							{
								rpo.getPrg()[curr][lact][mim][mpreg][sea] = rpo.getCow()[curr][lact][mim][mpreg][sea] - rpo.getCow()[curr][lact][mim][0][sea];
							//	System.out.println(" - rpo:%.0f = prg:%.0f",rpo.getcow[curr][lact][mim][0][sea],rpo.getprg[curr][lact][mim][mpreg][sea]);
							}
							else if(lact < bio.getMaxlact())	//late abortion, next lactation, dead calf not accounted for
							{
								rpo.getPrg()[curr][lact][mim][mpreg][sea] = rpo.getCow()[curr][lact][mim][mpreg][sea] - rpo.getCow()[curr][lact+1][1][0][sea];
							//	rpo.getprg[curr][lact][mim][mpreg][sea] = rpo.getcow[curr][lact][mim][mpreg][sea] - rpo.getcow[curr][lact][mim][0][sea];
							//	System.out.println(" - rpo:%.0f = prg:%.0f",rpo.getcow[curr][lact][mim][0][sea],rpo.getprg[curr][lact][mim][mpreg][sea]);
							}
							else if(lact == bio.getMaxlact())	//late abortion, no next lactation, dead calf not accounted for
							{
								rpo.getPrg()[curr][lact][mim][mpreg][sea] = rpo.getCow()[curr][lact][mim][mpreg][sea] - rpo.getCow()[curr][lact][mim][0][sea];
							//	rpo.getprg[curr][lact][mim][mpreg][sea] = rpo.getcow[curr][lact][mim][mpreg][sea] - rpo.getcow[curr][lact][mim][0][sea];
							//	System.out.println(" - rpo:%.0f = prg:%.0f",rpo.getcow[curr][lact][mim][0][sea],rpo.getprg[curr][lact][mim][mpreg][sea]);
							}
						}
					}
				}
			}
		}
		System.out.println("\n PREGNANCY_VALUE done ...");
		return;
	}/*** end PREGNANCY_VALUE ***/
	
	
	
	
	
}